# XTROEM-FV: A new code for computational astrophysics based on very high-order finite volume methods. Part I: 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. Part I: Magnetohydrodynamics. 2015. [hal-01199723v2](https://hal.science/hal-01199723v2)

HAL Id: hal-01199723

https://hal.science/hal-01199723v2

Preprint submitted on 27 Oct 2015

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_1.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.
## Part I: Magnetohydrodynamics

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

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

Tuesday 27th October, 2015

### ABSTRACT

The present work describes the building blocks of a new code for computational magnetohydrodynamics based on very high-order finite volume methods on Cartesian meshes. Spatial high-order accuracy is obtained with a weighted essentially non-oscillatory (WENO) reconstruction operator up to seventh order, while the time discretization is performed with a fourth order Strong-Stability Preserving Runge-Kutta method. Based on a shock detection approach, the reconstruction operator employs a very high-order WENO scheme in smooth flow regions and a third order WENO scheme in those parts of the flow with discontinuities or shocks. The Generalized Lagrange multiplier method is employed to enforce the solenoidal constraint on the magnetic field. Extensive numerical computations in one and two space dimensions are reported. Convergence rates for smooth flows verify the high-order accuracy of the scheme, and tests with strong shocks, including the Orszag-Tang vortex, the cylindrical blast wave problem, the rotor problem, and the Kelvin-Helmholtz instability, confirm the robustness and stability of the approach.

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

## 1 INTRODUCTION

The modeling of astrophysical phenomena has prompted the search for efficient and accurate numerical formulations for solving the magnetohydrodynamics (MHD) equations. Many codes have been developed for solving the MHD equations. Usually, the numerical methods chosen to be the main building block of these codes are the second order finite difference method (Ryu & Jones 1995a,b; Ryu et al. 1998), and the finite volume method (Zachary & Colella 1992; Zachary et al. 1994; Dai & Woodward 1994a,b; Balsara & Spicer 1999a; Tóth 2000; Janhunen 2000; Dedner et al. 2002; Ziegler 2004; Balsara 2004). Numerical frameworks for solving the magnetohydrodynamics equations with applications in astrophysics include the ZEUS CODE (Stone & Norman 1992a,b; Stone et al. 1992), the RAMSES CODE (Teyssier 2002), the PLUTO CODE (Mignone et al. 2007), the ATHENA CODE (Stone et al. 2008), and the ENZO CODE (Bryan et al. 2014). They make use of different reconstruction procedures, as well as adaptive mesh refinement algorithms in order to simulate accurately flow problems with complex and highly dynamical structures (Berger & Oliger 1984; Berger & Colella 1989; Bell et al. 1994; Balsara 2001; Tang & Tang 2003; Anderson et al. 2006; Zhang & MacFadyen 2006; Rosenberg et al. 2006; van der Holst

& Keppens 2007; Li 2010; Miniati & Martin 2011; Keppens et al. 2012; Hu 2013; Dumbser et al. 2013).

High-order numerical methods have gained quite popularity in the last years due to the need of high fidelity predictions in the simulations. Low-order methods show a considerable amount of numerical dissipation, meanwhile for high-order methods, this dissipation is very low. An advantage of low-order methods over high-order methods is their robustness in flows with discontinuities and shocks. Examples of high-order numerical methods for conservation laws are the conservative finite difference methods (Shu & Osher 1988, 1989), and the finite volume methods (Godunov 1959; van Leer 1979; Woodward & Colella 1984; Shu 2009). Both methods make use of a high-order reconstruction operator in order to achieve high-order accuracy in regions with smooth flows. The finite volume method has the advantage of working also on unstructured meshes, although the structure of the reconstruction operator is much more complicated as well as the selection of the stencil (Dumbser & Käser 2007; Dumbser et al. 2007). The conservative finite difference methods require uniform structured grids for the same purpose. The WENO schemes provide one of the most widely used reconstruction operators applied to solve numerically conservation laws. WENO stands for weighted essentially non-oscillatory and is based on the ENO schemes from Harten et al. (1987). In Liu et al. (1994) the WENO schemes were introduced, and in Jiang & Shu (1996) a general framework to construct arbitrary order accurate together with new smoothness indicators was provided. In Balsara & Shu

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

<footnote>† E-mail: munz@iag.uni-stuttgart.de</footnote>
(2000), smoothness indicators and interpolated values are reported up to seventh order of accuracy. In Gerolymos et al. (2009) very high-order WENO schemes are reported, up to 17th order of accuracy. In Dumbser & Käser (2007); Dumbser et al. (2007) are discussed WENO reconstruction operators for unstructured meshes for linear and nonlinear hyperbolic problems. In this work we construct a very high-order finite volume method on structured meshes for the magnetohydrodynamics equations that is robust and reliable for flows with very complex structures and for complicated shock configurations. The reconstruction is carried out with a WENO scheme.

Conservation laws allow in their solution discontinuities. It is well known that these discontinuities will produce spurious oscillations in the solution, the so-called Gibbs phenomenon (Gibbs 1898, 1899). These oscillations can cause unphysical states (negative density or pressure), making the code to fail. Therefore, a robust and efficient stabilization mechanism has to be constructed in order to get a well behaved numerical solution when shocks take place. For second-order finite volume and finite difference methods, this is mainly achieved by using slope/flux limiters (Harten et al. 1983; Harten 1983), and for high-order schemes, by using an adaptive reconstruction operator (like the WENO reconstruction (Jiang & Shu 1996; Shu 2009)). However, for some problems, the reconstruction can even produce these oscillations, hence it is advisable to perform the reconstruction on characteristic variables. We have observed that for systems of conservation laws, this last approach even fails for problems with strong shocks. For that reason, in this work we follow a simple, efficient and robust strategy based on the idea of shock detection, and further order reduction in the reconstruction step, but only applying this in the vicinity of the shock. This shock capturing approximation has been successfully employed in the past, and works very well, still for relativistic flows (Mignone et al. 2007; Tchekhovskoy et al. 2007; Beckwith & Stone 2011; Radice & Rezzolla 2011, 2012).

In the underlying base numerical scheme, the constraint $\nabla \cdot \mathbf{B} = 0$ is not satisfied. In fact, the spatial and time discretization produce errors in $\nabla \cdot \mathbf{B}$, which may grow with time, and may lead also to unphysical results (Brackbill & Barnes 1980; Powell 1994). For this reason, every MHD solver has to take into account a procedure to keep $\nabla \cdot \mathbf{B} = 0$. There are basically two approaches to overcome this problem: The divergence cleaning and the constrained transport algorithms. In the divergence cleaning approach, numerical errors in $\nabla \cdot \mathbf{B}$ are removed. The most popular schemes of this class are the Hodge projection (Chorin 1967; Brackbill & Barnes 1980), the 8-wave formulation (Powell 1994), and the Generalized Lagrange Multiplier (GLM) approach (Dedner et al. 2002). In the Hodge projection, the numerical solution of $\mathbf{B}$ is projected onto a subspace of zero divergence solutions. In spite of maintaining the solenoidal constraint up to machine accuracy, this scheme requires a time consuming solution of a Poisson equation at each time step. In the 8-wave formulation, an additional non-physical wave is introduced into the MHD system, giving as result the appearance of source terms proportional to the divergence of the magnetic field. In the Generalized Lagrange Multiplier approach, a new scalar field $\psi$ is introduced. This variable couples the divergence constraint with the evolution equation of the magnetic field. Local divergence errors are propagated to the boundary of the computational domain by a wave equation. This approach is very simple and its implementation is straightforward, and at the same time, conservation of all physical variables is maintained. From the other side, constrained transport algorithms exactly preserve the divergence constraint of the magnetic field, but with the price of introducing a staggered

mesh (Evans & Hawley 1988). The components of the magnetic field are defined at the cell interfaces and then interpolated to the cell barycenter. In this work we consider the GLM approach for controlling the solenoidal constraint of the magnetic field. In our algorithm we use the mixed hyperbolic/parabolic divergence cleaning as it is reported in Mignone et al. (2010).

The structure of this paper is as follows: We start by formulating the MHD equations in conservation form together with the GLM divergence cleaning for the preserving of the solenoidal constraint of the magnetic field. Next we describe in detail the multidimensional finite volume method with high-order WENO reconstruction and the shock capturing strategy used in this work. Finally, numerical computations of one-dimensional and two-dimensional problems with different shock configurations are shown and discussed.

## 2 IDEAL MAGNETOHYDRODYNAMICS

### 2.1 MHD Equations

Magnetohydrodynamics provides a powerful framework for describing the macroscopic behavior of plasmas, including both laboratory and space plasmas. The equations of the ideal magnetohydrodynamics, derived from the combination of Euler equations of the Hydrodynamics and the Maxwell equations of Electrodynamics, are given by the conservation of mass, the conservation of momentum, the conservation of energy and induction equations:

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

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

$$ \frac{\partial E}{\partial t} + \nabla \cdot (E\mathbf{v} + \mathbf{P} \cdot \mathbf{v}) = 0, \qquad (1c) $$

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

This system must satisfy an additional constraint: the solenoidal property of the magnetic field, which means that magnetic monopoles do not exist,

$$ \nabla \cdot \mathbf{B} = 0. \qquad (2) $$

The pressure tensor appearing in equations (1b) and (1c) combines the influence of the hydrodynamic and the magnetic pressure (the quantity in brackets is the total pressure)

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

An equation of state (EOS) is used to close the system. In this work we make use of the ideal gas equation of state with adiabatic exponent $\gamma$

$$ p = \rho \epsilon(\gamma - 1), \qquad (4) $$

where $\epsilon$ is the specific internal energy. For computational purposes, we rewrite the equation of state as a function of the total energy $E$, and the primitive variables, $\rho$, $v$, and $\mathbf{B}$

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

### 2.2 Divergence Cleaning with the Generalized Lagrange Multiplier Method

The divergence-free constraint of the magnetic field given in equation (2) has to be preserved also from the numerical point of view.
The simplest and straightforward manner to do it is by means of the divergence cleaning of Munz et al. (1999); Dedner et al. (2002). In this work we use the mixed hyperbolic/parabolic strategy first proposed in Dedner et al. (2002). In this approach, the solenoidal constraint (2) is coupled with the induction equation (1d) through the potential $\psi$. The induction equation and the solenoidal constraint take the following form

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

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

Equations (1a), (1b), (1c), (6a), and (6b), constitute the so-called GLM-MHD system. This equation system is still conservative (except for the equation for the unphysical scalar field $\psi$).

## 2.3 Eigenvalues of the MHD-GLM System

Next, we present the eigenvalues of the Jacobian matrices, $\mathbf{A}_\alpha$, for the GLM-MHD equations. These matrices are defined by

$$ \mathbf{A}_i = \frac{\partial \mathbf{f}_i(\mathbf{u})}{\partial \mathbf{u}}, \quad (i = x, y, z), \quad (7) $$

with $\mathbf{u}$ the state vector of conservative variables (plus the scalar field $\psi$) and $\mathbf{f} = [\mathbf{f}, \mathbf{g}, \mathbf{h}]$ the tensor of physical fluxes. The fluxes $\mathbf{f}, \mathbf{g}, \mathbf{h}$ are respectively the flux in the $x$-direction, in the $y$-direction, and in the $z$-direction. Here we show only the eigenvalues of the Jacobian matrices in the $x$-direction, $\mathbf{A}_x$, whereas the cases $y$ and $z$ easily follows from symmetry. The eigenvalues of $\mathbf{A}_x$ in nondecreasing order are

$$ \begin{aligned} \lambda_1 &= -c_h, & \lambda_2 &= v_x - c_f, & \lambda_3 &= v_x - c_a, \\ \lambda_4 &= v_x - c_s, & \lambda_5 &= v_x, & \lambda_6 &= v_x + c_s, \\ \lambda_7 &= v_x + c_a, & \lambda_8 &= v_x + c_f, & \lambda_9 &= c_h, \end{aligned} \quad (8) $$

where $c_f, c_a$, and $c_s$ are the fast, Alfvén, and slow characteristic speeds. Observe that in the equation (8) there are three MHD waves families associated with the characteristic speeds and an entropy mode. The three characteristic speeds are given by

$$ c_a = |b_1|, \quad (9) $$

$$ c_s = \sqrt{\frac{1}{2} \left( a^2 + b^2 - \sqrt{(a^2 + b^2)^2 - 4a^2 b_1^2} \right)}, \quad (10) $$

$$ c_f = \sqrt{\frac{1}{2} \left( a^2 + b^2 + \sqrt{(a^2 + b^2)^2 - 4a^2 b_1^2} \right)}, \quad (11) $$

where we have used the abbreviations

$$ a^2 = \frac{\gamma p}{\rho}, \quad b^2 = \frac{|\mathbf{B}|^2}{\rho}, \quad b_1^2 = \frac{B_x^2}{\rho}. \quad (12) $$

Observe the two new eigenvalues from the GLM-MHD equations, namely $\lambda_{1,9} = \pm c_h$. These eigenmodes propagate the divergence errors to the boundaries at speed $c_h$; besides, these divergence errors are damped at a rate $c_h^2/c_p^2$. The magnitude of the wave speeds $c_h$ is set to the maximum allowed speed in the pure MHD system and that is compatible with the CFL restriction in an explicit time discretization,

$$ c_h = \max_i (\max(|\lambda_{2,i}|, |\lambda_{8,i}|)), \quad i = x, y, z. \quad (13) $$

The value of the constant $c_p$ is chosen after setting $c_r \equiv c_p^2/c_h = 0.18$ (see Dedner et al. (2002)).

# 3 NUMERICAL METHODS

## 3.1 Finite Volume Methods

Let us consider the following system of conservation laws

$$ \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 (14) $$

By integrating equation (14) over the computational 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}}]$, we get the semi-discrete scheme

$$ \frac{\mathrm{d}\mathbf{u}_{ijk}}{\mathrm{d}t} = -\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}_{ij,k+\frac{1}{2}} - \hat{h}_{ij,k-\frac{1}{2}}}{\Delta z}, \quad (15) $$

where, in the context of finite volume methods, $\mathbf{u}_{ijk}$ is the spatial average of $\mathbf{u}$ in the cell $\Omega_{ijk}$ at time $t$

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

with $|\Omega_{ijk}| = \Delta x \Delta y \Delta z$, and $\hat{f}_{i \pm \frac{1}{2},jk}$, $\hat{g}_{i,j \pm \frac{1}{2},k}$, and $\hat{h}_{ij,k \pm \frac{1}{2}}$ are spatial averages of the physical fluxes over the cell faces $\bar{x}_{i \pm \frac{1}{2}}, y_{j \pm \frac{1}{2}}$, and $z_{k \pm \frac{1}{2}}$, respectively, at time $t$

$$ \begin{aligned} \hat{f}_{i \pm \frac{1}{2},jk} &= \frac{1}{|\sigma_{jk}|} \int_{\sigma_{jk}} \hat{f}(\mathbf{u}(x_{i \pm \frac{1}{2}}, y, z)) \, dz \, dy, \\ \hat{g}_{i,j \pm \frac{1}{2},k} &= \frac{1}{|\sigma_{ik}|} \int_{\sigma_{ik}} \hat{g}(\mathbf{u}(x, y_{j \pm \frac{1}{2}}, z)) \, dz \, dx, \end{aligned} \quad (17) $$

$$ \hat{h}_{ij,k \pm \frac{1}{2}} = \frac{1}{|\sigma_{ij}|} \int_{\sigma_{ij}} \hat{h}(\mathbf{u}(x, y, z_{k \pm \frac{1}{2}})) \, dy \, dx, $$

with the surface elements 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$. The numerical fluxes

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

are an approximation of the physical fluxes and are consistent 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); Toro (2009)).

When solving hyperbolic conservation laws, one of the most important requirement on the numerical method is that it should be a conservative scheme. In fact, the Lax-Wendroff theorem guarantees that for such numerical methods, if they are convergent, they converge to the weak solution (Lax & Wendroff 1960). Another important result related to conservative methods is that presented in (Hou & LeFloch 1994). They showed that a non-conservative numerical scheme does not converge to the correct solution if shock waves are present in the solution. Because of discontinuities are inherent to conservation laws, it is necessary to use conservative methods and those non-conservative should be avoided.

The semi-discrete scheme (15) is an exact relation as well as the averaged quantities (16) and (17). A high-order approximation of the fluxes (17) has to be found in order to get a high-order
4 Jonatan Núñez-de la Rosa, and Claus-Dieter Munz

accurate solution of (15). The integrals of the fluxes given in the equation (17) are discretized by means of a high-order Gaussian quadrature with suitable Gaussian integration points over the faces of the control volume; in the case of the present work, for Cartesian meshes, these are written as

$$
\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{19} \\
\hat{h}_{i j, 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}
$$

A high-order accurate numerical procedure is required to recon-
struct the point-wise values of **u** at the Gaussian integration points
at the faces from the only known quantities, the cell averages **u**<sub>ijk</sub> at
the cell barycenter. By evaluating the fluxes at the Gaussian points,
two sets of reconstructed values at a given face are present: those
obtained through the use of the cell Ω<sub>ijk</sub> as the main cell in the
reconstruction procedure and those associated with the neighbor-
ing cell. These are known as the left and right values at the face
interface: **u**<sub>L</sub> and **u**<sub>R</sub>. Finally the fluxes are evaluated by replacing
**f**, **g**, and **h** by a monotone flux, the so-called Riemann solver (Toro
2009).

XTROEM-FV includes several Riemann solvers for the MHD equations, namely, the Rusanov numerical flux (Rusanov 1962), the Roe Riemann solver (Roe & Balsara 1996), the HLL Riemann solver (Toro 2009), the HLLC Riemann solver (Toro et al. 1994; Gurski 2004), and the HLLD Riemann solver (Miyoshi & Kusano 2005). In spite of the better accuracy of the HLL family of Riemann solvers, in this work we report computation that only use the Rusanov Riemann solver because it is much more robust and cheap that the other ones (because it adds enough dissipation to ensure stability of the scheme (LeVeque 1992)). The Rusanov numerical flux is given by

$$
\hat{f}(u_L, u_R) = \frac{1}{2} (f(u_L) + f(u_L) - |\lambda_{\max}|(u_L - u_R)), \quad (20)
$$

where $\lambda_{\max}$ is the largest local wave speed, which guarantees the
stability of the scheme. Taking $\lambda_{\max}$ locally instead of globally
implies that less viscosity will be added in regions where the solution
is smooth.

Finally, a flow diagram for the XTROEM-FV framework based
on finite volume methods is depicted in the figure 1.

3.2 Discretization of the GLM Divergence Cleaning

The equations for $B_x$ and $\psi$ are decoupled from the rest of the GLM-MHD system. In fact, cleaning a vector quantity $\mathbf{B}$ with divergence errors results in 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}. \qquad (21)
$$

The numerical flux is then derived as the solution of the local
Riemann problem with left-hand state $(B_{x,l}, \psi_l)^T$ and right hand
state $(B_{x,r}, \psi_r)^T$ as

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

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

<figure><img src="image_2.png" /><figcaption>Figure 1. Flow diagram for the high-order finite volume method.</figcaption></figure>

It is possible to employ for the quantities $B_x$ and $\psi$ the Riemann
solver used in the full MHD system, but in this work we use the
exact solution (22) of the linear Riemann problem (21). Dedner et al.
(2002) suggest to use the solution (22) as input for the Riemann
solver used in the solution of the other conserved quantities.

There are basically two possible ways to deal with the source
term in equation (6b). The first one consists in simple add this source
term to the hyperbolic update, that is to the semi-discrete scheme
(15) (see Susanto et al. (2013)). The second one is based on an
operator-splitting approach. Following the idea presented in Dedner
et al. (2002), we first solve the homogeneous GLM-MHD system in
a so-called *hyperbolic step*, and then we consider the source term
in the *source step*. The scalar field *ψ* 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 (24)
$$

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

The basic idea of the WENO schemes is based on an adaptive
reconstruction procedure to obtain a higher-order approximation
on smooth regions while the scheme remains non-oscillatory near
discontinuities (Shu 2009). In each stencil, a polynomial is recon-
structed from the cell averages of the solution and later a weighted combination of all these polynomials is constructed. The nonlinear weights are computed taking into account the smoothness of every polynomial in its respective stencil.

In the following, the dimension-by-dimension WENO reconstruction algorithm is outlined. In appendix A is presented the description of the main aspects of the WENO reconstruction operator, as for example, the computation of the oscillation indicators, and the linear and nonlinear weights for the WENO3, WENO5 and WENO7 schemes. The dimension-by-dimension reconstruction was introduced in Casper & Atkins (1993) in the ENO framework for two-dimensional conservation laws. The scheme was extended to three-dimensional domains in Titarev & Toro (2004). We emphasize that the dimension-by-dimension reconstruction is only valid for structured meshes, and at the same time, it is the less computationally expensive among all reconstructions (although a new methodology for finite volume methods in Cartesian meshes has been introduced in Buchmüller & Helzel (2014), making the reconstruction process as simple as for finite difference methods). For unstructured grids (and they are also valid in Cartesian grids) there are another strategies, for example, the genuine multidimensional reconstruction of Hu & Shu (1999), but they will not be considered in this work.

The aim of the reconstruction is the following: From the cell averages of the solution, points values of the function at appropriate points (that is, the Gaussian integration points) are interpolated (or reconstructed). These will be used in the finite-volume methodology to calculate the numerical fluxes via a Riemann solver. This means, given the cell averages<fnref n="1" /> of the function $u(x, y, z)$ in the cell $\Omega_{ijk}$

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

we reconstruct 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}})$. In this work we make use of the two points Gaussian quadrature rule

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

for calculating the surface integral appearing in the finite volume formulation just as it is reported in Titarev & Toro (2004, 2005). In the figure 2 these Gaussian integration points are depicted. They 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}$. A higher-order Gaussian quadrature gives as result negatives linear weights in the point-wise WENO reconstruction. Although there is an strategy to deal with such negatives weights (see Shi et al. (2002)), we prefer to follow the recommendations given in Titarev & Toro (2004). We stress that the high-order accuracy of the scheme is provided by the high-order reconstruction of the function values at the Gaussian integration points. In spite of this two-points quadrature leads to a formal fourth order of accuracy, several numerical experiments confirm that the order of accuracy is not damaged, and also the computations are more efficient (Titarev & Toro 2004).

In the following, we are going to show the dimension-by-dimension algorithm only for the reconstruction at the Gaussian

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Gaussian integration points in a control volume where the WENO reconstruction is performed. Left: At cell edges in the two-dimensional case, they are given by $(x_{i\pm\frac{1}{2}}, y_{j+\alpha})$ and $(x_{i+\alpha}, y_{j\pm\frac{1}{2}})$. Right: At cell faces in the three-dimensional case, they are given by $(x_{i\pm\frac{1}{2}}, y_{j+\alpha}, z_{k+\beta})$, $(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}})$. In both cases $\alpha, \beta = \pm 1/2\sqrt{3}$.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Stencil for the high-order WENO5 reconstruction for finite volume methods. The state is reconstructed in the Gaussian integration points depicted in the Figure 2. Left: For the cell $\Omega_{ij}$ in 2D. Right: For the cell $\Omega_{ijk}$ in 3D.</figcaption></figure>

integration points $(x_{i\pm\frac{1}{2}}, y_{j+\alpha}, z_{k+\beta})$, i.e., the points at faces $x_{i\pm\frac{1}{2}}$. In an analog way we can reconstruct 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}})$, corresponding to the faces $y_{j\pm\frac{1}{2}}$ and $z_{k\pm\frac{1}{2}}$. The algorithm consists of three sweeps in the 3D case and two sweeps in the 2D case. Let us start defining the stencils (and adopting the notation given in Titarev & Toro (2004)). In order 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 stencil are of degree $N$), we require the stencil to be formed by the cells $\Omega_{i_x, i_y, i_z}$, with $i_x, i_y$, and $i_z$ satisfying

$$\begin{align}
i - N \le i_x \le i + N, \nonumber \\
j - N \le i_y \le j + N, \nonumber \\
k - N \le i_z \le k + N. \tag{27}
\end{align}$$

As an example, the stencils for the WENO5 reconstruction procedure are depicted in the Figure 3 for the two- and three-dimensional cases. Once the stencils are defined, we proceed with the description of the sweeps:

(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 indexes $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

<footnote n="1">¹ The quantity $u$ is a dummy function, and it can be interpreted as a primitive or a conservative variable. Its cell average is given by $u_{ijk}$. It should not be confused $u$ as a dummy function, and $u$, the vector of conservative variables.</footnote>
6 Jonatan Núñez-de la Rosa, and Claus-Dieter Munz

<figure><img src="image_5.png" /><figcaption>Figure 4. Sweeps and stencils for the two-dimensional WENO5 reconstruction procedure. Top: x-direction. Bottom: y-direction.</figcaption></figure>

$\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} \begin{vmatrix} x_{i \pm \frac{1}{2}} \\ y_{j \pm \frac{1}{2\sqrt{3}}} \end{vmatrix} = \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.

The dimension-by-dimension WENO reconstruction for the two-dimensional case can be summarized as follows (see figure 4)

$$ \{u_{ij}\} \rightarrow \{\bar{u}_{i \pm \frac{1}{2}, i_y}\} \rightarrow \{u_{i \pm \frac{1}{2}, j \pm \frac{1}{2\sqrt{3}}}\}, \quad (28) $$

and for the three-dimensional case

$$ \{u_{ijk}\} \rightarrow \{\bar{u}_{i \pm \frac{1}{2}, i_y} \}_{i_z} \rightarrow \{\bar{u}_{i \pm \frac{1}{2}, j \pm \frac{1}{2\sqrt{3}}, i_z}\} \quad (29) \\ \rightarrow \{u_{i \pm \frac{1}{2}, j \pm \frac{1}{2\sqrt{3}}, k \pm \frac{1}{2\sqrt{3}}}\}. $$

In XTROEM-FV, the WENO reconstruction operator of order third, fifth, and seventh order are implemented. In this work, the reported two-dimensional computations were done with the WENO7 operator.

## 3.4 Shock Capturing for High-Order Finite Volume Methods

The general methodology used in this work to stabilize the solution when the scheme fails to give a physically meaningful state, because either shocks are present or some unphysical states are produced in the intermediate stages of the simulation (e.g., during the reconstruction step), consists basically in the three blocks outlined in the following

(1) Checking of quantities with restrictions

* Positivity of density

* Positivity of pressure

(2) Detection of regions with strong shocks

* Shocks indicators

* Marking of troubled points/cells and direct neighbors

(3) Special treatment in troubled regions

* Employment of robust second/third order scheme

* Employment of robust Riemann solvers

The way this strategy is implemented establishes different “flavors” of the shock capturing, but in principle, the ground idea is the shock detection and a further order reduction of the scheme in the vicinity of the shock, the so-called *fallback approaches* (Mignone et al. 2007; Tchekhovskoy et al. 2007; Beckwith & Stone 2011; Radice & Rezzolla 2011, 2012). This is equivalent to adding more dissipation, because a lower-order scheme together with a dissipative Riemann solver is employed in the vicinity of the shock. There are some cases in which the high-order WENO reconstruction fails when it is performed, that is, negative densities or pressure are obtained after the reconstruction step. This occurs especially in regions with severe shocks or when in the stencil several discontinuities are present. In such cases, the WENO algorithm cannot select an oscillation-free polynomial from all stencils. For that reason, we have to adopt the reconstruction-order-reduction strategy outlined above in those stencils containing more than one discontinuity or having a strong shock. That is, we first apply a shock indicator to every point/cell and then mark only those having a shock and their direct neighbors. For this purpose, we use 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}|}. \quad (30) $$

If the Jameson indicator is larger than $\eta = 5.0 \times 10^{-3}$, then the cell is flagged. Next, for every marked cell, a WENO3 reconstruction (or in severe cases, the second order MUSCL reconstruction) is used. If these schemes fail to give a physical state, then we switch to the first-order Godunov scheme, where no reconstruction is performed. We remark that the use of the above strategy will never deteriorate the high order of accuracy of the method for smooth solutions, and the reason is that the indicator only flags the cells with discontinuities. We refer to the classical work of Harten et al. (1987) for more details.

We have to point out that this approach has provided very good results for a wide series of problems in one-dimensional and multidimensional magnetohydrodynamics. In fact, all calculations presented in this work have made use of WENO3 reconstruction in troubled cells and higher-order WENO schemes in smooth parts of the flow.

## 3.5 Time Discretization

The semi-discrete scheme (15) is solved explicitly using the so-called *method of lines*. The method of lines acquires the accuracy order of the integrator used to solve the system of ordinary differential equations, under the condition that the spatial discretization is of the same order of accuracy or higher (Schiesser 1991). The system of ordinary differential equations (15) is solved with the family of Strong-Stability Preserving Runge-Kutta Methods (SSPRK).
The class of SSPRK methods was first developed in Shu & Osher (1988); Shu (1988), where they called them TVD (Total Variation Diminishing) time discretizations, and further by Gottlieb & Shu (1998). New developments have been reported in Hundsdorfer et al. (2003); Ruuth & Spiteri (2002); Spiteri & Ruuth (2002, 2003); Gottlieb (2005); Gottlieb et al. (2009). The main idea of these methods resides in assuming that the first-order forward Euler method is strongly stable under the total variation norm (and in general, any given norm) and a suitable time step restriction. From this, the aim is to construct a higher-order time discretization that preserves strong stability under the TVD norm, and maybe, with a different time step restriction (Shu & Osher 1988). The SSPRK methods are desirable in problems with discontinuities and strong shocks, because they guarantee that, as part of the time-integration process, no additional oscillations are introduced (Hesthaven & Warburton 2008; Shu 2009).

The systems of ordinary differential equations (15) can be recast as

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

where $-L(u, t)$ is the spatial discretization operator. Let us consider the general m-stages Runge-Kutta method of the Shu-Osher form (Shu & Osher 1988)

$$ \begin{aligned} u^{(0)} &= u^n, \\ u^{(i)} &= \sum_{k=0}^{i-1} \left( \alpha_{ik} u^{(k)} + \Delta t \beta_{ik} L(u^{(k)}, t^n + \gamma_k \Delta t) \right), \quad (32) \\ u^{n+1} &= u^{(m)}, \end{aligned} $$

with $\alpha_{ik} \ge 0$, $i = 1, \dots, m$. The coefficients of the fourth-order SSPRK with five stages derived in Spiteri & Ruuth (2002) are presented in the table 1. This scheme has an SSP coefficient $c = \min_{i,k} \frac{\alpha_{ik}}{\beta_{ik}} = 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.

### 3.5.1 Courant-Friedrichs-Lewy Condition

The Courant-Friedrichs-Lewy (CFL) condition is a necessary condition for stability of any explicit one-level numerical scheme (Courant et al. 1928). It establishes that the domain of dependence of the solution is contained in the numerical domain of dependence of the numerical method. From a physical point of view, this condition guarantees that the propagation speed of any traveling wave is always smaller than the numerical speed $c_N \equiv \Delta x/\Delta t$. In other words, the time step should be equal to or smaller than the time necessary for a wave to travel across the stencil of the spatial discretization scheme (Blazek 2005). This condition is applied to constrain the time step, and in the one-dimensional case is given by

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

with $\lambda_k$ the maximum eigenvalue of the physical x-flux Jacobian over all computational cells. The number $C_{CFL}$ satisfies $C_{CFL} \le 1$, and depends on the spatial discretization and on the time-stepping scheme used for solving the conservation law. In the multidimensional case, only approximated values of the CFL condition are available for nonlinear equations. We consider here the CFL condition presented in Titarev & Toro (2005), namely

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

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), \qquad (35) $$

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

## 4 NUMERICAL COMPUTATIONS

### 4.1 Propagation of a smooth circularly polarized Alfvén wave

This problem is characterized for being an exact solution of the MHD equations. Because of the smoothness of this wave, it can be used for doing a convergence analysis. The test is run with the following simulation parameters (Tóth 2000): The dimensions of the computational domain satisfy $L_x = 2L_y = 1.0$, and the discretization of it, $N_x = 2N_y$. The wave propagates along the diagonal of the grid, at an angle $\alpha = \arctan(0.5) \approx 26.6^\circ$ with respect to the x-axis. Since the wave does not propagate along the diagonals of the grid cells, it is guaranteed that the x- and y-fluxes are different; as a consequence, the problem is truly multi-dimensional. The density is set to $\rho = 1$, and the pressure $p = 0.1$. Let be $v_0 = 0.1$ and $B_0 = 0.1$. Additionally let be $x_\parallel = (x \cos \alpha + y \sin \alpha)$. The components of the velocity and magnetic field perpendicular to the wave vector are given by $v_\perp = v_0 \sin(2\pi x_\parallel)$ and $B_\perp = B_0 \sin(2\pi x_\parallel)$, and the velocity and magnetic field in z-direction are $v_z = v_0 \cos(2\pi x_\parallel)$ and $B_z = B_0 \cos(2\pi x_\parallel)$. The quantities $v_\perp$ and $B_\perp$ are related to $B_x$ and $B_y$ via $B_\perp = B_y \cos \alpha - B_x \sin \alpha$, and $B_\parallel = B_x \cos \alpha + B_y \sin \alpha$. An adiabatic equation of state with $\gamma = 5/3$ is considered. At the boundaries of the domain, periodic boundary conditions are imposed. In the table 2 are shown the convergence rates at time $t_f = 1.0$ when the MHD equations are solved with the finite volume scheme using the WENO reconstruction operator. The experimental order of convergence (EOC) confirms the high-order accuracy of the numerical method for smooth flows, even when nonlinearities are present in the exact solution.

### 4.2 One-dimensional Riemann Problems

In this section we present the calculation of one-dimensional Riemann problems, which feature different shocks configurations after the initial state is evolved to certain end time $t_f$. The study of shock tube problems constitutes the standard benchmark in computational magnetohydrodynamics because it is through these simple initial conditions that the robustness of the scheme is evaluated and also it is analyzed its capacity to resolve continuous and discontinuous flow structures (Mignone & Tzeferacos 2010). The considered tests are those discussed in Ryu & Jones (1995a), although some of them were first presented in Brio & Wu (1988) and in Dai & Woodward (1994b). For all these tests, we set the adiabatic index $\gamma = 5/3$ and transmissive boundary conditions. The CFL number is set to $C_{CFL} = 0.95$. The one-dimensional domain is the closed interval $[0, 1]$, and the membrane separating the left and right initial states is localized in the point $x = 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$. Observe that left and right states are provided in primitive variables. The computational domain is decomposed into 500 cells. The WENO3, WENO5, and WENO7 reconstruction operators were used in smooth parts of the flow and the WENO3 in regions with discontinuities and/or shocks (only for the WENO5 and WENO7 schemes). The six initial conditions are given in table 3. The exact solutions have been obtained
<table><caption>**Table 1.** Coefficients of the fourth-order SSPRK method with five stages.</caption><thead><tr><th>Stages</th><th></th><th>k = 0</th><th>k = 1</th><th>k = 2</th><th>k = 3</th><th>k = 4</th></tr></thead><tbody><tr><td rowspan="5">α<sub>ik</sub></td><td>i = 1</td><td>1</td><td></td><td></td><td></td><td></td></tr><tr><td>i = 2</td><td>0.44437049406734</td><td>0.55562950593266</td><td></td><td></td><td></td></tr><tr><td>i = 3</td><td>0.62010185138540</td><td>0</td><td>0.37989814861460</td><td></td><td></td></tr><tr><td>i = 4</td><td>0.17807995410773</td><td>0</td><td>0</td><td>0.82192004589227</td><td></td></tr><tr><td>i = 5</td><td>0.00683325884039</td><td>0</td><td>0.51723167208978</td><td>0.12759831133288</td><td>0.34833675773694</td></tr><tr><td rowspan="5">β<sub>ik</sub></td><td>i = 1</td><td>0.39175222700392</td><td></td><td></td><td></td><td></td></tr><tr><td>i = 2</td><td>0</td><td>0.36841059262959</td><td></td><td></td><td></td></tr><tr><td>i = 3</td><td>0</td><td>0</td><td>0.25189177424738</td><td></td><td></td></tr><tr><td>i = 4</td><td>0</td><td>0</td><td>0</td><td>0.54497475021237</td><td></td></tr><tr><td>i = 5</td><td>0</td><td>0</td><td>0</td><td>0.08460416338212</td><td>0.22600748319395</td></tr><tr><td>γ<sub>k</sub></td><td></td><td>0</td><td>0.39175222700392</td><td>0.58607968896779</td><td>0.47454236302687</td><td>0.93501063100924</td></tr></tbody></table>

<table><caption>**Table 2.** Convergence rates for the MHD equations with initial condition given by the smooth circularly polarized Alfvén wave. In these calculations we have employed the Rusanov Riemann solver. Results for $L_2$ norm error of the y-component of magnetic field are given, with reconstruction operator WENO3, WENO5 and WENO7. Simulation time was set to $t_f = 1.0$.</caption><thead><tr><th>Method</th><th>Cells</th><th>L<sub>2</sub> error</th><th>L<sub>2</sub> order</th></tr></thead><tbody><tr><td rowspan="5">WENO3</td><td>20 × 10</td><td>2.163 × 10<sup>-3</sup></td><td></td></tr><tr><td>40 × 20</td><td>3.218 × 10<sup>-4</sup></td><td>2.75</td></tr><tr><td>80 × 40</td><td>4.120 × 10<sup>-5</sup></td><td>2.97</td></tr><tr><td>160 × 80</td><td>5.097 × 10<sup>-6</sup></td><td>3.01</td></tr><tr><td>320 × 160</td><td>6.275 × 10<sup>-7</sup></td><td>3.02</td></tr><tr><td rowspan="5">WENO5</td><td>20 × 10</td><td>6.375 × 10<sup>-5</sup></td><td></td></tr><tr><td>40 × 20</td><td>2.134 × 10<sup>-6</sup></td><td>4.90</td></tr><tr><td>80 × 40</td><td>6.848 × 10<sup>-8</sup></td><td>4.96</td></tr><tr><td>160 × 80</td><td>2.158 × 10<sup>-9</sup></td><td>4.98</td></tr><tr><td>320 × 160</td><td>6.561 × 10<sup>-11</sup></td><td>5.03</td></tr><tr><td rowspan="5">WENO7</td><td>20 × 10</td><td>8.482 × 10<sup>-6</sup></td><td></td></tr><tr><td>40 × 20</td><td>8.285 × 10<sup>-8</sup></td><td>6.67</td></tr><tr><td>80 × 40</td><td>7.984 × 10<sup>-10</sup></td><td>6.69</td></tr><tr><td>160 × 80</td><td>7.543 × 10<sup>-12</sup></td><td>6.72</td></tr><tr><td>320 × 160</td><td>6.308 × 10<sup>-14</sup></td><td>6.90</td></tr></tbody></table>

with the exact Riemann solver developed by Torrilhon (2002, 2003), and they are plotted as red lines. The description of the flow patterns is presented below together with some comments regarding the capabilities of the XTROEM-FV CODE for solving such problems.

**Riemann Problem 1 (RP-1)** This initial state will develop into two fast shocks, one weak slow shock, one slow rarefaction and one contact discontinuity. Plots of the density, thermal pressure, and the y-component of the magnetic field are depicted in the figure 5 (top). The final simulation time is $t_f = 0.08$. Small oscillations are observed only for the WENO7 scheme. The shocks are very well captured in around 3-4 cells with the WENO3 reconstruction operator. The WENO5 scheme provides the best results in this simulation.

**Riemann Problem 2 (RP-2)** The final simulation time is $t_f = 0.03$. This Riemann problem features one left-moving fast shock and one left-moving slow shock, one fast rarefaction and one slow rarefaction wave moving to the right, and one contact discontinuity. The z-component of the velocity and magnetic field is zero in all

regions. Plots of the density, thermal pressure, and the y-component of the magnetic field are depicted in the figure 5 (middle). Some oscillations are present for the WENO5 and WENO7 schemes. In spite of that, the shocks are better resolved when they are combined with the WENO3 scheme.

**Riemann Problem 3 (RP-3)** This problem is characterized by the appearance of multiple weak discontinuities, namely, two fast shocks, two slow shocks, two rotational discontinuities and one contact discontinuity. The shocks propagate from each side of the contact discontinuity. Plots of the density, thermal pressure, and the y-component of the magnetic field at time $t_f = 0.20$ are depicted in the figure 5 (bottom). The WENO5 and WENO7 solutions required less cells (around 3-4) to resolve the discontinuities than the WENO3 scheme (around 10-20).

**Riemann Problem 4 (RP-4)** This test involves one of the so-called switch-on and switch-off structures, a switch-on fast shock. The main property of this problems is that the tangential magnetic field turns on in the region located behind the switch-on fast shocks and switch-on slow rarefactions waves. The switch-on fast shock is moving to the right. A fast rarefaction, a slow rarefaction, a contact discontinuity, and a slow shock are among the other structures developed in this problem. Plots of the density, thermal pressure, and the y-component of the magnetic field are shown in the figure 6 (top) at the time $t_f = 0.15$. Only the solution computed with WENO7 scheme shows oscillations around the discontinuities, which are resolved using only 3-4 cells.

**Riemann Problem 5 (RP-5)** This problem comprises a switch-on slow rarefaction wave. In this case, the following structures are created: A hydrodynamic rarefaction, a switch-on slow rarefaction, a contact discontinuity, a slow shock, a rotational discontinuity, and a fast rarefaction. Plots of the density, thermal pressure, and the y-component of the magnetic field at time $t_f = 0.16$ are depicted in the figure 6 (middle). The structures mentioned before are showed in the figure, from left to right. The shocks have been very well captured using around 3-4 cells with the WENO3 reconstruction.

**Riemann Problem 6 (RP-6)** This problem is similar to the Sod problem for the Euler equations. The test is characterized by the formation of the so-called compound structures. These structures involve a shock and a rarefaction wave traveling together. This kind
<table><caption>**Table 3.** Riemann problems for the magnetohydrodynamics equations. The initial states have been taken from Ryu & Jones (1995a). The adiabatic index is set to $\gamma = 5/3$ and the computational domain is the interval $[0, 1]$. For the problems 1 to 3, the magnetic field components have to be divided by $\sqrt{4\pi}$.</caption><thead><tr><th>Test</th><th>State</th><th>ρ</th><th>ν<sub>x</sub></th><th>ν<sub>y</sub></th><th>ν<sub>z</sub></th><th>p</th><th>B<sub>x</sub></th><th>B<sub>y</sub></th><th>B<sub>z</sub></th></tr></thead><tbody><tr><td rowspan="2">RP-1</td><td>w<sub>L</sub></td><td>1.00</td><td>10.00</td><td>0.00</td><td>0.00</td><td>20.00</td><td>5.00</td><td>5.00</td><td>0.00</td></tr><tr><td>w<sub>R</sub></td><td>1.00</td><td>-10.00</td><td>0.00</td><td>0.00</td><td>1.00</td><td>5.00</td><td>5.00</td><td>0.00</td></tr><tr><td rowspan="2">RP-2</td><td>w<sub>L</sub></td><td>1.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>1.00</td><td>3.00</td><td>5.00</td><td>0.00</td></tr><tr><td>w<sub>R</sub></td><td>0.10</td><td>0.00</td><td>0.00</td><td>0.00</td><td>10.00</td><td>3.00</td><td>2.00</td><td>0.00</td></tr><tr><td rowspan="2">RP-3</td><td>w<sub>L</sub></td><td>1.08</td><td>1.20</td><td>0.01</td><td>0.50</td><td>0.95</td><td>2.00</td><td>3.60</td><td>2.00</td></tr><tr><td>w<sub>R</sub></td><td>1.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>1.00</td><td>2.00</td><td>4.00</td><td>2.00</td></tr><tr><td rowspan="2">RP-4</td><td>w<sub>L</sub></td><td>1.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>1.00</td><td>1.00</td><td>1.00</td><td>0.00</td></tr><tr><td>w<sub>R</sub></td><td>0.20</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.10</td><td>1.00</td><td>0.00</td><td>0.00</td></tr><tr><td rowspan="2">RP-5</td><td>w<sub>L</sub></td><td>1.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>1.00</td><td>0.70</td><td>0.00</td><td>0.00</td></tr><tr><td>w<sub>R</sub></td><td>0.30</td><td>0.00</td><td>0.00</td><td>1.00</td><td>0.20</td><td>0.70</td><td>1.00</td><td>0.00</td></tr><tr><td rowspan="2">RP-6</td><td>w<sub>L</sub></td><td>1.00</td><td>0.00</td><td>0.00</td><td>0.00</td><td>1.00</td><td>0.75</td><td>1.00</td><td>0.00</td></tr><tr><td>w<sub>R</sub></td><td>0.125</td><td>0.00</td><td>0.00</td><td>0.00</td><td>0.10</td><td>0.75</td><td>-1.00</td><td>0.00</td></tr></tbody></table>

of structures was first analyzed in Brio & Wu (1988). This test involves a left-going slow compound wave. Plots of the density, thermal pressure, and the y-component of the magnetic field at time $t_f = 0.1$ are shown in the figure 6 (bottom). The plots show the development (from left to right) of a left-going fast rarefaction, the left-going slow compound wave, a contact discontinuity, a right-going slow shock wave, and a right-going fast rarefaction wave. The WENO3 reconstruction operator used for shock capturing resolves very well the shock wave and the contact discontinuity within just a few cells.

## 4.3 Multidimensional Test Problems

### 4.3.1 Magnetic Field Loop Advection

The magnetic field loop advection is a very important test for multidimensional MHD because with it, 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. Because of the loop remains in magnetostatic balance, after some periods, its profile should be the same. With multidimensional MHD codes that do not include an algorithm for the treatment of the solenoidal constraint of the magnetic field, or if they use numerical methods too much diffusive, the magnetic loop will smear over the time. More details can be found in Tóth (2000); Stone et al. (2008); Mignone & Tzeferacos (2010).

Following the description to this problem given in Tóth (2000), the computational domain is given by 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 domain. The initial velocity is given by

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

where $v_0$ is the magnitude of the velocity (we use $v_0 = \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$ (we use $R = 0.3$). In this way, for $r \le R$ we have

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

where $r = \sqrt{x^2 + y^2}$, $B_0$ is the magnitude of the magnetic field and it should be small in such a way that the magnetic pressure is smaller than the gas pressure (we use $B_0 = 10^{-3}$). An adiabatic equation of

state with $\gamma = 5/3$ is considered. At the boundaries of the domain we consider periodic boundary conditions. The simulation time is $t_f = 2.0$, and in the figure 7 is depicted the magnetic pressure at times $t \in \{0.0, 2.0\}$ when the RKFV method was used. The computational domain is decomposed into $800 \times 400$ cells. The WENO7 reconstruction is used. No shock capturing method was employed in 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. The CFL condition was set to $C_{CFL} = 0.95$. We observe that the loop profile is only kept when the GLM is taken into account.

### 4.3.2 Current Sheet

This problem was first discussed in Hawley & Stone (1995). Initially, a region is uniformly filled with a gas at rest. The magnetic field is initialized in such a way that this switches signs at the slices $x = +0.25$ and $x = -0.25$. Then we perturb the system with a sinusoidal velocity function in y, which generates nonlinear, linearly polarized Alfvén waves. Because the magnetic pressure does not remain constant, these Alfvén waves turn into magnetosonic waves. Magnetic reconnection occurs because of the two current sheets at $x = \pm 0.25$. Additionally, since the parameter $\beta < 1$, the magnetic reconnection drives highly over-pressurized regions, which launch magnetosonic waves transverse to the field, causing magnetic energy to be transformed into thermal energy (Hawley & Stone 1995). Large magnetic field gradients are produced close to the points where the magnetic reconnection take place. Because of that, this test problem is useful to check if the algorithm can handle these features.

For the simulation, we consider the computational domain to be the box $[-0.5, +0.5] \times [-0.5, +0.5]$. We set the density and pressure uniform in the whole domain, with $\rho = 1.0$, and $p = 0.5\beta$, where $\beta$ is an input parameter which represents the ratio of gas pressure to magnetic energy density. We set $v_y = v_z = 0$, and $v_x = A \sin(2\pi y)$, where $A$ is a parameter that is typically used to test the robustness of the algorithm. We set $B_x = B_z = 0$, and $B_y = 1$ for $|x| > 0.25$ and $B_y = -1$ otherwise. We employ $A = 0.1$ and $\beta = 0.1$ in our simulations. An adiabatic equation of state with $\gamma = 5/3$ is considered. At the boundaries of the domain we consider periodic boundary conditions. The simulation time is $t_f = 10.0$, and in the figure 8 are depicted the density at times $t \in \{5.0, 7.5\}$ obtained with the finite volume scheme with WENO7 reconstruction. The CFL condition was set to $C_{CFL} = 0.95$. The results are quite similar
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Plots for the Riemann problems RP-1, RP-2, and RP-3. The initial states have been provided in the table 3. An ideal equation of state with adiabatic index $\gamma = 5/3$ is used. The computations were done with a finite volume method with WENO3, WENO5, and WENO7 reconstruction operators. The CFL number is $C_{CFL} = 0.95$.</figcaption></figure>

to those obtained in Stone et al. (2008), but also they show the robustness of the code when tackling complex flow configurations with very high-order WENO reconstruction operators.

### 4.3.3 Orszag-Tang Vortex

This problem was first studied in Orszag & Tang (1979) for the incompressible MHD equations. Many authors have used the Orszag-Tang vortex for the compressible MHD equations as a test problem in order to know how robust is the employed numerical scheme at handling the formation and the interactions of MHD shocks (Zachary et al. 1994; Ryu & Jones 1995a; 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, which consists of smooth initial data, is obtained by the superposition of a velocity vortex with a magnetic vortex. Because of this highly unstable setup, a

broad range of MHD waves are generated, interacting with each other, making a transition towards turbulence.

For this problem, the computational domain is given by the box $[0, 1] \times [0, 1]$. We set the density and pressure uniform in the whole domain, with $\rho = \gamma^2$ and $p = \gamma$. With this choice of the density and pressure we have a sound speed $c_s = 1$. For the 2D case of the Orszag-Tang vortex, the velocity has the following profile

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

The magnetic field is given by

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

An adiabatic equation of state with $\gamma = 5/3$ is considered. At the boundaries of the domain we consider periodic boundary conditions. The simulation time is $t_f = 1.0$, and in the figure 9 are depicted the density and magnetic pressure at time $t = 0.5$. The computational domain is decomposed into 600 × 600 cells. The
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Solution plots for the Riemann problems RP-4, RP-5, and RP-6. The initial states have been provided in the table 3. An ideal equation of state with adiabatic index $\gamma = 5/3$ is used. The computations were done with a finite volume method employing the WENO3, WENO5, and WENO7 reconstruction operators. The CFL number is $C_{CFL} = 0.95$.</figcaption></figure>

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 used, with $c_r = 0.18$, and $c_h$ determined by the maximum propagation speed in the system. We set the CFL condition to $C_{CFL} = 0.95$. For this very demanding problem, all schemes (WENO3, WENO5 and WENO7) are very stable and robust. The code was able to handle very well the shock waves and their interactions with the other flow structures emerging in the evolution of this configuration. In the figure 10 is plotted the pressure along the slices $y = 0.4277$ and $y = 0.3125$ at time $t = 0.5$ for different mesh resolutions and different WENO reconstruction operators. It is observed the lack of oscillations in these slices for all WENO reconstructions.

#### 4.3.4 Cylindrical Blast Wave

This problem concerns with the formation and propagation of strong shock waves in a highly magnetized medium (Zachary et al. 1994; Londrillo & del Zanna 2000; Stone et al. 2008; Mignone et al. 2010). Initially, a cylindrical region located in the center of a domain is filled with a magnetized overpressured gas. As a result, a strong shock wave moving outwards is formed. This test is very difficult for any numerical method in the sense that the code can crash because of the formation of unphysical values in quantities like the density or the pressure. If periodic boundary conditions are used, then the interactions of the shock waves lead to very complex configurations.

For this problem, the computational domain is given by the box $[0, 1] \times [0, 1]$. We set the density and pressure uniform in the whole domain, with $\rho = 1$ and $p = 0.1$. 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 = 10.0$. Initially, the velocity is set to zero, that is $v_x = v_y = v_z = 0$. For the 2D case, the magnetic field
<figure><img src="image_8.png" /><figcaption>**Figure 7.** 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 = 2.0$ (right). The solution was computed with a finite volume method with WENO7 reconstruction. The computational domain is the box $[-1.0, +1.0] \times [-0.5, +0.5]$, discretized by using a mesh of 800 × 400 cells.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Current sheet problem. Contour plots of the density at times $t = 5.0$ (left) and $t = 7.5$ (right). The calculations were performed with a RKFV 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>

is set to $B_x = B_y = B_z = 1/\sqrt{2}$. An adiabatic equation of state with $\gamma = 5/3$ is considered. At the boundaries of the domain we consider periodic boundary conditions. The simulation time is $t_f = 0.5$, and in the figure 11 are depicted the density and magnetic pressure at times $t = 0.1$, and $t = 0.5$. The computational domain is decomposed into 600 × 600 cells. The reconstruction operator used for this problem was the WENO7, combined with the WENO3 for shock regions. 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$. The finite volume scheme effectively resolves the strong gradients arising after the release of the cylindrical region.

### 4.3.5 Rotor Problem

The rotor problem was first proposed in Balsara & Spicer (1999b), and also it was studied in Tóth (2000). It consists of a high-density, rapidly spinning fluid in a low-density fluid. Initially, both fluids are subject to a uniform magnetic field. Because of the rapidly rotating fluid, torsional Alfvén waves are launched into the fluid at rest. As a consequence the rotor suffers a decrease in its angular momentum.

The initial configuration of this problem is as follows: The computational domain is the box $[0, 1] \times [0, 1]$. We set the density and pressure uniform in the ambient medium, with $\rho = 1.0$ and $p = 1.0$. Inside 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)$) we set the density $\rho = 10.0$, and the
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Orszag-Tang vortex problem. Contour plots of the density (left) and the magnetic pressure (right) at time $t = 0.5$. Computations were done with a RKFV method 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_11.png" /><figcaption>**Figure 10.** Orszag-Tang vortex problem. Plots of the pressure along the slices $y = 0.3125$ (top) and $y = 0.4277$ (bottom) at time $t = 0.5$ for different mesh resolutions and different WENO reconstruction operators.</figcaption></figure>

pressure is the same as in the ambient fluid. The velocity in the ambient medium is initially set to zero, that is $v_x = v_y = v_z = 0$. The rotor has an angular velocity $\omega$ such that $v = \omega r = 1$ at $r = 0.1$. A linear taper is applied to the velocity and 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 set to $B_x = 5.0$, and $B_y = B_z = 0$ in the whole computational domain. An adiabatic equation of state with $\gamma = 7/5$ is considered. At the boundaries of the domain we consider transmissive boundary conditions. The computational domain is discretized using 600 × 600 finite volume cells. The 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, and $c_h$ was determined by the maximum wave propagation speed in the system. We set the CFL condition to $C_{CFL} = 0.95$.

This problem is extremely difficult for any numerical method, even for low-order schemes. A mechanism to preserve the positivity of the density and pressure has to be kept in mind for this problem, and efficiently combined with the shock capturing approach. The end simulation time is $t_f = 0.5$, and in the figure 12 are depicted the density, pressure, Mach number, and the magnitude of the magnetic field at time $t = 0.25$. In the plot for the Mach number is clearly visible that the fluid inside the rotor is still rotating with uniform angular speed. In the figure 13 is plotted the pressure along the slices $x = 0.5$ and $y = 0.5$ at time $t = 0.25$ for different mesh resolutions and different WENO reconstruction operators.

### 4.3.6 Cloud-Shock Interaction

This problem consists in the interaction of a strong shock wave with a dense cloud. We follow mainly the setup presented by Tóth (2000), which is basically similar to the discussed by Dai & Woodward (1994b).

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 given by

$$ w_L = \begin{pmatrix} 3.86859 \\ 0.0 \\ 0.0 \\ 0.0 \\ 167.345 \\ 0.0 \\ 2.1826182 \\ -2.1826182 \end{pmatrix}, \quad w_R = \begin{pmatrix} 1.0 \\ -11.2536 \\ 0.0 \\ 0.0 \\ 1.0 \\ 0.0 \\ 0.56418958 \\ 0.56418958 \end{pmatrix}. \qquad (40) $$

The cloud is assumed to be a high-density cylinder, which center is located at $(1.6, 0.5)$. Its radius is $r = 0.15$, and its density and pressure are $\rho = 10.0$, $p = 1.0$. The cloud is in hydrostatic equilibrium with the ambient gas. An adiabatic equation of state with $\gamma = 5/3$ is considered. At the boundaries of the domain we consider transmissive boundary conditions. The end simulation time is $t_f = 0.1$,
<figure><img src="image_12.png" /><figcaption>Figure 11. Cylindrical blast wave. Contour plots of the density (left) and magnetic pressure (right) at times $t = 0.10$ (top) and $t = 0.50$ (bottom) obtained with a finite volume method 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>

and in the figure 14 are depicted the density, gas pressure, and the magnitude of the velocity at time $t = 0.1$, when the MHD equations were solved with the RKFV method and using the WENO7 reconstruction with WENO3 as shock capturing. The CFL condition was set to $C_{CFL} = 0.95$.

### 4.3.7 Kelvin-Helmholtz Instability

The Kelvin-Helmholtz instability occurs when a perturbation is introduced 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 nonlinear magnetohydrodynamic turbulence. As a test of the linear regime, one can compare the growth rate of the instability with the analytic result before the instability becomes nonlinear. A single mode perturbation is needed for such a comparison.

The computational domain is the box $[-0.5, +0.5] \times [-0.5, +0.5]$. The density is set to $\rho = 1.0$ for $|y| \le 0.25$, and $\rho = 2.0$ for $|y| < 0.25$. We set the pressure uniform in the whole

domain, with $p = 2.5$. We set $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 (41)$$

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 (42)$$

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 equation of state with $\gamma = 7/5$ is considered. At the boundaries of the domain we consider periodic boundary conditions. The computational domain is decomposed into 600 × 600 cells. The WENO7 reconstruction is used in smooth parts of the flow and the WENO3 in regions with discontinuities. The divergence cleaning of Dedner et al. (2002) was used. The CFL condition was set to $C_{CFL} = 0.95$. The simulation is run until $t_f = 6.0$, and in the fig-
<figure><img src="image_13.png" /><figcaption>**Figure 12.** Rotor problem. Contour plots of the density, pressure, Mach number, and the magnitude of the magnetic field at time $t = 0.25$. These calculations were performed with a 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>

ure 15 are depicted the density at times $t \in \{1.0, 2.0, 3.0, 6.0\}$. In the figure 16 is plotted the area-averaged of the y-component of the velocity during the linear growth phase of the Kelvin-Helmholtz instability when a high-order finite volume scheme with different resolutions and WENO reconstruction operators is used.

# 5 CONCLUSIONS

XTROEM-FV is a new simulation code for computational astrophysics based on state-of-the-art very high-order finite volume WENO schemes. In this work, we have focused on the implementation details of the ideal magnetohydrodynamics equations in XTROEM-FV. Algorithms for achieving arbitrary high-order of accuracy in space are discussed in detail on Cartesian meshes. For this purpose, the WENO schemes are used in combination with a fallback approach as shock capturing strategy. The point values of the solution at cell interfaces are obtained with the reconstruction procedure. The input values of the reconstruction operator are the cell averages

of the solution written in primitive variables. It is not necessary to use the characteristic variables because the use of the fallback approach, and as consequence, the expensive computation of the left and right eigenvectors of the physical flux Jacobian matrices is avoided. This scheme is quite robust and offers the possibility of using arbitrary high-order reconstruction operators in those smooth parts of the flow, and stable reconstruction operators (typically second or third order accurate) when discontinuities or strong rarefactions take place.

XTROEM-FV includes several Riemann solvers for computing the solution values at cell interfaces. The Rusanov numerical flux showed to be the most robust and stable among all numerical fluxes considered in our calculations. The other Riemann solvers implemented (HLL, HLLC and HLLD) failed to compute the cell interface values for some test problems with very high-order of accuracy, especially those involving strong shock waves and complicated interactions of these with the other flow structures. Strong Stability-Preserving Runge-Kutta schemes were also used in order to integrate in time the system of equations obtained from the finite volume dis-
<figure><img src="image_14.png" /><figcaption>**Figure 13.** Rotor problem. Plots of the y-component of the magnetic field along the slice $y = 0.5$ (top) and plots of the x-component of the magnetic field along the slice $x = 0.5$ (bottom) at time $t = 0.25$ for different mesh resolutions and different WENO reconstruction operators.</figcaption></figure>

cretization. The divergence constraint of the magnetic field in the MHD equations has to be kept from the numerical point of view. For that purpose, the divergence cleaning mechanism of Dedner et al. (2002) has been used because its simplicity, ease of implementation, and robustness. High-order numerical computations confirm the reliability and feasibility of the divergence cleaning for keeping the solenoidal constraint.

XTROEM-FV was subject to a wide range of problems in one and two space dimensions. Multidimensional smooth flow problems were used for probing the high-order of accuracy of the numerical algorithms, and the results confirm that XTROEM-FV numerically converges to the expected rate. Additionally, one-dimensional and two-dimensional flows containing very challenging structures like discontinuities, strong rarefactions, and shock waves, have been simulated with XTROEM-FV. The followed strategy for the shock capturing allowed us to compute such problems with very high-order WENO reconstruction operators. All two-dimensional calculations reported in this work were obtained with a WENO7 operator, combined with a WENO3 for troubled zones. Computations with WENO5 and WENO7 reconstruction operators look very similar for this kind of problems, although only those obtained with WENO7 have been reported.

Further work involves the implementation on XTROEM-FV of the equations of relativistic astrophysics on flat spacetimes (relativistic hydrodynamics and relativistic magnetohydrodynamics), resistive magnetohydrodynamics and resistive relativistic magnetohydrodynamics. Moreover, numerical flow simulations with complex and highly dynamical structures require efficient adaptive mesh refinement (AMR) algorithms (see for example, Berger & Oliger (1984); Berger & Colella (1989); Bell et al. (1994); Balsara (2001); Teyssier (2002); Anderson et al. (2006); Zhang & MacFadyen (2006); van der Holst & Keppens (2007); Miniati & Martin (2011); Keppens et al. (2012); Hu (2013); Dumbser et al. (2013); Bryan et al. (2014)), because such methods can save a lot of computational effort and memory consumption over a static mesh approach. A new version of

<figure><img src="image_15.png" /><figcaption>**Figure 14.** Cloud-Shock Interaction. Contour plots of the density (top), pressure (middle), and velocity magnitude (bottom) at time $t = 0.1$. The computations were done with a RKFV with WENO7 reconstruction. The computational domain is the box $[0, 2] \times [0, 1]$, discretized by using a mesh of 800 $\times$ 400 cells.</figcaption></figure>

XTROEM-FV including AMR techniques is under development. We are focusing on unstructured hexahedral meshes with tree-based adaptive mesh refinement algorithms, which allow recursive grid refinements on a cell-by-cell basis (Teyssier 2002; van der Holst & Keppens 2007; Keppens et al. 2012). With AMR, XTROEM-FV will be able to control the mesh resolution and resolve the small flow structures in an very efficient manner. Applications to magnetohydrodynamic turbulence will be also considered within this context.

## ACKNOWLEDGMENTS

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 and Partho Upadhyay for the very careful reading of the manuscript and their valuable suggestions, which improved the quality.
<figure><img src="image_16.png" /><figcaption>**Figure 15.** Kelvin-Helmholtz instability. From left to right and from top to bottom, contour plots of the density at simulation times $t \in \{1.0, 2.0, 3.0, 6.0\}$. The solution was computed with a RKFV-WENO7 scheme on a mesh made of 600 $\times$ 600 cells.</figcaption></figure>

## REFERENCES

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

* Balsara D., 2001, J. Comput. Phys., 174, 614

* Balsara D., 2004, Astrophys. J. Suppl. Ser., 151, 149

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

* Balsara D., Spicer D., 1999a, J. Comput. Phys., 148, 133

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

* Beckwith K., Stone J., 2011, Astrophys. J. Suppl. Ser., 193, 6

* Bell J., Berger M., Saltzman J., Welcome M., 1994, SIAM J. Sci. Comput., 15, 127

* Berger M., Colella P., 1989, J. Comput. Phys., 82, 64

* Berger M., Oliger J., 1984, J. Comput. Phys., 53, 484

* 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

* Bryan G., et al., 2014, Astrophys. J. Suppl. Ser., 211, 19

* Buchmüller P., Helzel C., 2014, J. Sci. Comput., pp 1–26

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

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

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

* Dai W., Woodward P., 1994a, J. Comput. Phys., 111, 354

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

* Dai W., Woodward P., 1998, Astrophys. J., 494, 317

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

* 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., Zanotti O., Hidalgo A., Balsara D., 2013, J. Comput. Phys., 248, 257

* Evans C., Hawley J., 1988, Astrophys. J., 332, 659

* Gerolymos G., Sénéchal D., Vallet I., 2009, J. Comput. Phys., 228, 8481

* Gibbs J. W., 1898, Nature, 59, 200

* Gibbs J. W., 1899, Nature, 59, 606

* 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
<figure><img src="image_17.png" /><figcaption>**Figure 16.** Kelvin-Helmholtz instability. Area-averaged of the y-component of the velocity during the linear growth phase of the Kelvin-Helmholtz instability. Computations were done with a high-order finite volume scheme with different resolutions and WENO reconstruction operators.</figcaption></figure>

Gurski K., 2004, SIAM J. Sci. Comput., 25, 2165

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

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-Verlag, Heidelberg

Hou T., LeFloch P., 1994, Math. Comp., 62, 497

Hu G., 2013, J. Comput. Phys., 252, 591

Hu C., Shu C.-W., 1999, J. Comput. Phys., 150, 97

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

Janhunen P., 2000, J. Comput. Phys., 160, 649

Jiang G.-S., Shu C.-W., 1996, J. Comput. Phys., 126, 202

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

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

Li S., 2010, J. Comput. Appl. Math., 233, 3139

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

Londrillo P., del Zanna L., 2000, Astrophys. J., 530, 508

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

Mignone A., Bodo G., Massaglia S., Matsakos T., Tesileanu O., Zanni C., Ferrari A., 2007, Astrophys. J. Suppl. Ser., 170, 228

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

Miniati F., Martin D., 2011, Astrophys. J. Suppl. Ser., 195, 5

Miyoshi T., Kusano K., 2005, J. Comput. Phys., 208, 315

Munz C.-D., Schneider R., Sonnendrücker E., Voss U., 1999, C. R. Acad. Sci. Math., 328, 431

Orszag S., Tang C.-M., 1979, J. Fluid Mech., 90, 129

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, Astron. Astrophys., 547, A26

Roe P., Balsara D., 1996, SIAM J. Appl. Math., 56, 57

Rosenberg D., Fournier A., Fischer P., Pouquet A., 2006, J. Comput. Phys., 215, 59

Rusanov V., 1962, USSR Comput. Math. Phys., 1, 304

Ruuth S., Spiteri R., 2002, J. Sci. Comput., 17, 211

Ryu D., Jones T., 1995a, Astrophys. J., 442, 228

Ryu D., Jones T., 1995b, Astrophys. J., 452, 785

Ryu D., Miniati F., Jones T., Frank A., 1998, Astrophys. J., 509, 244

Schiesser W., 1991, The numerical method of lines: Integration of partial differential equations. Academic Press

Shi J., Hu C., Shu C.-W., 2002, J. Comput. Phys., 175, 108

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., Norman M., 1992a, Astrophys. J. Suppl. Ser., 80, 753

Stone J., Norman M., 1992b, Astrophys. J. Suppl. Ser., 80, 791

Stone J., Mihalas D., Norman M., 1992, Astrophys. J. Suppl. Ser., 80, 819

Stone J., Gardiner T., Teuben P., Hawley J., Simon J., 2008, Astrophys. J. Suppl. Ser., 178, 137

Susanto A., Ivan L., Sterck H. D., Groth C., 2013, J. Comput. Phys., 250, 141

Tang H., Tang T., 2003, SIAM J. Numer. Anal., 41, 487

Tchekhovskoy A., McKinney J., Narayan R., 2007, Mon. Not. R. Astron. Soc., 379, 469

Teyssier R., 2002, Astron. Astrophys., 385, 337

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-Verlag, Berlin

Toro E., Spruce M., Speares W., 1994, Shock Waves, 4, 25

Torrilhon M., 2002, Technical report, Exact solver and uniqueness conditions for Riemann problems of ideal magnetohydrodynamics. Seminar for Applied Mathematics, ETH

Torrilhon M., 2003, J. Plasma Phys., 69, 253

Tóth G., 2000, J. Comput. Phys., 161, 605

Woodward P., Colella P., 1984, J. Comput. Phys., 54, 115

Zachary A., Colella P., 1992, J. Comput. Phys., 99, 341

Zachary A., Malagoli A., Colella P., 1994, SIAM J. Sci. Comput., 15, 263

Zhang W., MacFadyen A., 2006, Astrophys. J. Suppl. Ser., 164, 255

Ziegler U., 2004, J. Comput. Phys., 196, 393

van Leer B., 1979, J. Comput. Phys., 32, 101

van der Holst B., Keppens R., 2007, J. Comput. Phys., 226, 925

# APPENDIX A: WENO RECONSTRUCTION PROCEDURE

## A1 Basic Idea

We are going to present the WENO reconstruction of third, fifth and seventh order for the following interpolation points:

$$x = x_{i-\frac{1}{2}}, \quad x = x_{i+\frac{1}{2}}, \quad x = x_{i-\frac{1}{2\sqrt{3}}}, \quad x = x_{i+\frac{1}{2\sqrt{3}}}. \qquad (A1)$$

The general idea will be described in the following for the interpolation point $x = x_{i+\frac{1}{2}}$ but the procedure is completely similar for the other points.

Let $I = [a, b]$. Let $a < \dots < x_{i-1} < x_i < x_{i+1} < \dots < b$ be a partition of $I$. Let us assume that the mesh is uniform, i.e., $\Delta x = x_{i+1} - x_i = \text{constant}$. Therefore we can take $x_i = i\Delta x$. Let $u(x)$ be a function. Let us assume that all cell averages

$$\bar{u}_i = \frac{1}{\Delta x} \int_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}} u(x) \, dx \qquad (A2)$$
over the intervals $I_i = (x_{i-\frac{1}{2}}, x_{i+\frac{1}{2}})$ are given. As an example, let us find an approximation of the function $u(x)$ at a point other than the cell barycenter $x_i$, for instance, the point $x = x_{i+\frac{1}{2}}$. Defining the primitive function $u(x)$ by

$$U(x) = \int_{x-\frac{1}{2}}^{x} u(\zeta) \, d\zeta \qquad (A3)$$

where the lower limit is irrelevant and can be replaced for any fixed point, then clearly we have

$$\begin{aligned} U(x_{i+\frac{1}{2}}) &= \int_{x-\frac{1}{2}}^{x_{i+\frac{1}{2}}} u(\zeta) \, d\zeta \\ &= \sum_{l=0}^{i} \int_{x_{l-\frac{1}{2}}}^{x_{l+\frac{1}{2}}} u(\zeta) \, d\zeta = \sum_{l=0}^{i} \Delta x \bar{u}_l \end{aligned} \qquad (A4)$$

That is, knowing all the cell averages $\bar{u}_l$, we also know the point values of the primitive function $U(x_{i+\frac{1}{2}})$ at all half nodes. Therefore, interpolation polynomials can be constructed for the primitive function $U(x)$. The derivative of such an interpolation polynomial for $U(x)$ can then be used as an approximation to $u(x) = U'(x)$.

Let $P_k(x)$ be the polynomials of degree at most $N+1$ which interpolates the function $U(x)$ at the $N+2$ points $x_{j+\frac{1}{2}}$, for $j = i - N - 3 + k + l$, with $l = 1, \dots, N+2$. Let us approximate the value of the function $u(x)$ at $x_{i+\frac{1}{2}}$, $u_{i+\frac{1}{2}}$ by using polynomials of degree $N$. For that, we can use $N+1$ possible stencils $S_k$ containing the interval $I_i$ and build unique interpolation polynomials of degree at most $N$, $p_k(x) = P_k'(x)$ that reconstruct the function $u(x)$ over the stencil $S_k = \{I_j | j = i - N - 2 + k + l, \text{ with } l = 1, \dots, N+1\}$ in the sense that

$$(\bar{p}_k)_j = \frac{1}{\Delta x} \int_{x_{j-\frac{1}{2}}}^{x_{j+\frac{1}{2}}} p_k(x) \, dx = \bar{u}_j, \qquad (A5)$$

approximates the value $u_{i+\frac{1}{2}}^{(k)} \equiv p_k(x_{i+\frac{1}{2}})$, for $j = i - N - 2 + k + l$, with $l = 1, \dots, N+1$. All these polynomials satisfy $p_k(x_i) = u_i$.

## A2 WENO3 scheme

Let us approximate the value of the function $u(x)$ at the points $x_{i-\frac{1}{2}}$, $x_{i+\frac{1}{2}}$, $x_{i-\frac{1}{2\sqrt{3}}}$ and $x_{i+\frac{1}{2\sqrt{3}}}$ by using polynomials of degree $N=1$. For the stencil $S_k$ we get the corresponding values $u_{i-\frac{1}{2}}^{(k)}$, $u_{i+\frac{1}{2}}^{(k)}$, $u_{i-\frac{1}{2\sqrt{3}}}^{(k)}$ and $u_{i+\frac{1}{2\sqrt{3}}}^{(k)}$. Let us write down the stencils

$$S_1 = \{I_{i-1}, I_i\}, \quad S_2 = \{I_i, I_{i+1}\}. \qquad (A6)$$

The corresponding interpolated values for the point $x_{i-\frac{1}{2}}$ are

$$u_{i-\frac{1}{2}}^{(1)} = \frac{1}{2}(\bar{u}_{i-1} + \bar{u}_i), \qquad (A7)$$

$$u_{i-\frac{1}{2}}^{(2)} = \frac{1}{2}(3\bar{u}_i - \bar{u}_{i+1}), \qquad (A8)$$

and the corresponding interpolated values for the point $x_{i+\frac{1}{2}}$ are

$$u_{i+\frac{1}{2}}^{(1)} = \frac{1}{2}(-\bar{u}_{i-1} + 3\bar{u}_i), \qquad (A9)$$

$$u_{i+\frac{1}{2}}^{(2)} = \frac{1}{2}(\bar{u}_i + \bar{u}_{i+1}). \qquad (A10)$$

In the same way as before, the corresponding interpolated values for the point $x_{i-\frac{1}{2\sqrt{3}}}$ are

$$u_{i-\frac{1}{2\sqrt{3}}}^{(1)} = \frac{1}{6}(\sqrt{3}\bar{u}_{i-1} + 6\bar{u}_i - \sqrt{3}\bar{u}_i), \qquad (A11)$$

$$u_{i-\frac{1}{2\sqrt{3}}}^{(2)} = \frac{1}{6}(\sqrt{3}\bar{u}_i + 6\bar{u}_i - \sqrt{3}\bar{u}_{i+1}), \qquad (A12)$$

and the corresponding interpolated values for the point $x_{i+\frac{1}{2\sqrt{3}}}$ are

$$u_{i+\frac{1}{2\sqrt{3}}}^{(1)} = \frac{1}{6}(-\sqrt{3}\bar{u}_{i-1} + 6\bar{u}_i + \sqrt{3}\bar{u}_i), \qquad (A13)$$

$$u_{i+\frac{1}{2\sqrt{3}}}^{(2)} = \frac{1}{6}(-\sqrt{3}\bar{u}_i + 6\bar{u}_i + \sqrt{3}\bar{u}_{i+1}). \qquad (A14)$$

If we choose the large stencil

$$S = \{I_{i-1}, I_i, I_{i+1}\} \qquad (A15)$$

which is the union of all 2 stencils $S_k$, then we are able to find an interpolation polynomial $p(x)$ of degree at most 2, satisfying $p(x_j) = u_j$ for $i - N \le j \le i + N$ and giving the approximations $u_{i-\frac{1}{2}} \equiv p(x_{i-\frac{1}{2}})$, $u_{i+\frac{1}{2}} \equiv p(x_{i+\frac{1}{2}})$,

$$u_{i-\frac{1}{2}} = \frac{1}{6}(2\bar{u}_{i-1} + 5\bar{u}_i - \bar{u}_{i+1}), \qquad (A16)$$

$$u_{i+\frac{1}{2}} = \frac{1}{6}(-\bar{u}_{i-1} + 5\bar{u}_i + 2\bar{u}_{i+1}), \qquad (A17)$$

and the approximations $u_{i-\frac{1}{2\sqrt{3}}} \equiv p(x_{i-\frac{1}{2\sqrt{3}}})$ and $u_{i+\frac{1}{2\sqrt{3}}} \equiv p(x_{i+\frac{1}{2\sqrt{3}}})$

$$u_{i-\frac{1}{2\sqrt{3}}} = \frac{1}{12}(\sqrt{3}\bar{u}_{i-1} + 12\bar{u}_i - \sqrt{3}\bar{u}_{i+1}), \qquad (A18)$$

$$u_{i+\frac{1}{2\sqrt{3}}} = \frac{1}{12}(-\sqrt{3}\bar{u}_{i-1} + 12\bar{u}_i + \sqrt{3}\bar{u}_{i+1}), \qquad (A19)$$

provided that the function is smooth in the large stencil $S$.

## A3 WENO5 scheme

Let us approximate the value of the function $u(x)$ at the points $x_{i-\frac{1}{2}}$, $x_{i+\frac{1}{2}}$, $x_{i-\frac{1}{2\sqrt{3}}}$ and $x_{i+\frac{1}{2\sqrt{3}}}$ by using polynomials of degree $N=2$. For the stencil $S_k$ we get the corresponding values $u_{i-\frac{1}{2}}^{(k)}$, $u_{i+\frac{1}{2}}^{(k)}$, $u_{i-\frac{1}{2\sqrt{3}}}^{(k)}$ and $u_{i+\frac{1}{2\sqrt{3}}}^{(k)}$. Let us write down the stencils

$$\begin{aligned} S_1 &= \{I_{i-2}, I_{i-1}, I_i\}, \\ S_2 &= \{I_{i-1}, I_i, I_{i+1}\}, \end{aligned} \qquad (A20)$$

$$S_3 = \{I_i, I_{i+1}, I_{i+2}\}. \qquad (A21)$$

The corresponding interpolated values for the point $x_{i-\frac{1}{2}}$ are

$$u_{i-\frac{1}{2}}^{(1)} = \frac{1}{6}(-\bar{u}_{i-2} + 5\bar{u}_{i-1} + 2\bar{u}_i), \qquad (A21)$$

$$u_{i-\frac{1}{2}}^{(2)} = \frac{1}{6}(2\bar{u}_{i-1} + 5\bar{u}_i - \bar{u}_{i+1}), \qquad (A22)$$

$$u_{i-\frac{1}{2}}^{(3)} = \frac{1}{6}(11\bar{u}_i - 7\bar{u}_{i+1} + 2\bar{u}_{i+2}), \qquad (A23)$$
and the corresponding interpolated values for the point $x_{i+\frac{1}{2}}$ are

$$u_{i+\frac{1}{2}}^{(1)} = \frac{1}{6}(2\bar{u}_{i-2} - 7\bar{u}_{i-1} + 11\bar{u}_{i}), \qquad (A24)$$

$$u_{i+\frac{1}{2}}^{(2)} = \frac{1}{6}(-\bar{u}_{i-1} + 5\bar{u}_{i} + 2\bar{u}_{i+1}), \qquad (A25)$$

$$u_{i+\frac{1}{2}}^{(3)} = \frac{1}{6}(2\bar{u}_{i} + 5\bar{u}_{i+1} - \bar{u}_{i+2}). \qquad (A26)$$

In the same way as before, the corresponding interpolated values for the point $x_{i-\frac{1}{2\sqrt{3}}}$ are

$$u_{i-\frac{1}{2\sqrt{3}}}^{(1)} = \frac{1}{12}(-\sqrt{3}\bar{u}_{i-2} + 4\sqrt{3}\bar{u}_{i-1} + 12\bar{u}_{i} - 3\sqrt{3}\bar{u}_{i}), \qquad (A27)$$

$$u_{i-\frac{1}{2\sqrt{3}}}^{(2)} = \frac{1}{12}(\sqrt{3}\bar{u}_{i-1} + 12\bar{u}_{i} - \sqrt{3}\bar{u}_{i+1}), \qquad (A28)$$

$$u_{i-\frac{1}{2\sqrt{3}}}^{(3)} = \frac{1}{12}(12\bar{u}_{i} + 3\sqrt{3}\bar{u}_{i} - 4\sqrt{3}\bar{u}_{i+1} + \sqrt{3}\bar{u}_{i+2}), \qquad (A29)$$

and the corresponding interpolated values for the point $x_{i+\frac{1}{2\sqrt{3}}}$ are

$$u_{i+\frac{1}{2\sqrt{3}}}^{(1)} = \frac{1}{12}(\sqrt{3}\bar{u}_{i-2} - 4\sqrt{3}\bar{u}_{i-1} + 12\bar{u}_{i} + 3\sqrt{3}\bar{u}_{i}), \qquad (A30)$$

$$u_{i+\frac{1}{2\sqrt{3}}}^{(2)} = \frac{1}{12}(-\sqrt{3}\bar{u}_{i-1} + 12\bar{u}_{i} + \sqrt{3}\bar{u}_{i+1}), \qquad (A31)$$

$$u_{i+\frac{1}{2\sqrt{3}}}^{(3)} = \frac{1}{12}(12\bar{u}_{i} - 3\sqrt{3}\bar{u}_{i} + 4\sqrt{3}\bar{u}_{i+1} - \sqrt{3}\bar{u}_{i+2}). \qquad (A32)$$

If we choose the large stencil

$$S = \{I_{i-2}, I_{i-1}, I_i, I_{i+1}, I_{i+2}\} \qquad (A33)$$

which is the union of all 3 stencils $S_k$, then we are able to find an interpolation polynomial $p(x)$ of degree at most 4, satisfying $p(x_j) = u_j$ for $i-N \le j \le i+N$ and giving the approximations $u_{i-\frac{1}{2}} \equiv p(x_{i-\frac{1}{2}})$, $u_{i+\frac{1}{2}} \equiv p(x_{i+\frac{1}{2}})$,

$$u_{i-\frac{1}{2}} = \frac{60}{12}(-3\bar{u}_{i-2} + 27\bar{u}_{i-1} + 47\bar{u}_{i} - 13\bar{u}_{i+1} + 2\bar{u}_{i+2}), \qquad (A34)$$

$$u_{i+\frac{1}{2}} = \frac{1}{60}(+2\bar{u}_{i-2} - 13\bar{u}_{i-1} + 47\bar{u}_{i} + 27\bar{u}_{i+1} - 3\bar{u}_{i+2}), \qquad (A35)$$

and the approximations $u_{i-\frac{1}{2\sqrt{3}}} \equiv p(x_{i-\frac{1}{2\sqrt{3}}})$ and $u_{i+\frac{1}{2\sqrt{3}}} \equiv p(x_{i+\frac{1}{2\sqrt{3}}})$

$$u_{i-\frac{1}{2\sqrt{3}}} = \frac{1}{4320}(-70\sqrt{3}\bar{u}_{i-2} - \bar{u}_{i-2} + 500\sqrt{3}\bar{u}_{i-1} + 4\bar{u}_{i-1} + 4314\bar{u}_{i} - 500\sqrt{3}\bar{u}_{i+1} + 4\bar{u}_{i+1} + 70\sqrt{3}\bar{u}_{i+2} - \bar{u}_{i+2}), \qquad (A36)$$

$$u_{i+\frac{1}{2\sqrt{3}}} = \frac{1}{4320}(+70\sqrt{3}\bar{u}_{i-2} - \bar{u}_{i-2} - 500\sqrt{3}\bar{u}_{i-1} + 4\bar{u}_{i-1} + 4314\bar{u}_{i} + 500\sqrt{3}\bar{u}_{i+1} + 4\bar{u}_{i+1} - 70\sqrt{3}\bar{u}_{i+2} - \bar{u}_{i+2}), \qquad (A37)$$

provided that the function is smooth in the large stencil S.

## A4 WENO7 scheme

Let us approximate the value of the function $u(x)$ at the points $x_{i-\frac{1}{2}}$, $x_{i+\frac{1}{2}}$, $x_{i-\frac{1}{2\sqrt{3}}}$ and $x_{i+\frac{1}{2\sqrt{3}}}$ by using polynomials of degree $N=3$. For

the stencil $S_k$ we get the corresponding values $u_{i-\frac{1}{2}}^{(k)}$, $u_{i+\frac{1}{2}}^{(k)}$, $u_{i-\frac{1}{2\sqrt{3}}}^{(k)}$ and $u_{i+\frac{1}{2\sqrt{3}}}^{(k)}$. Let us write down the stencils

$$S_1 = \{I_{i-3}, I_{i-2}, I_{i-1}, I_i\}, \quad S_2 = \{I_{i-2}, I_{i-1}, I_i, I_{i+1}\}, \\ S_3 = \{I_{i-1}, I_i, I_{i+1}, I_{i+2}\}, \quad S_4 = \{I_i, I_{i+1}, I_{i+2}, I_{i+3}\}. \qquad (A38)$$

The corresponding interpolated values for the point $x_{i-\frac{1}{2}}$ are

$$u_{i-\frac{1}{2}}^{(1)} = \frac{1}{12}(\bar{u}_{i-3} - 5\bar{u}_{i-2} + 13\bar{u}_{i-1} + 3\bar{u}_{i}), \qquad (A39)$$

$$u_{i-\frac{1}{2}}^{(2)} = \frac{1}{12}(-\bar{u}_{i-2} + 7\bar{u}_{i-1} + 7\bar{u}_{i} - \bar{u}_{i+1}), \qquad (A40)$$

$$u_{i-\frac{1}{2}}^{(3)} = \frac{1}{12}(3\bar{u}_{i-1} + 13\bar{u}_{i} - 5\bar{u}_{i+1} + \bar{u}_{i+2}), \qquad (A41)$$

$$u_{i-\frac{1}{2}}^{(4)} = \frac{1}{12}(25\bar{u}_{i} - 23\bar{u}_{i+1} + 13\bar{u}_{i+2} - 3\bar{u}_{i+3}), \qquad (A42)$$

and the corresponding interpolated values for the point $x_{i+\frac{1}{2}}$ are

$$u_{i+\frac{1}{2}}^{(1)} = \frac{1}{12}(-3\bar{u}_{i-3} + 13\bar{u}_{i-2} - 23\bar{u}_{i-1} + 25\bar{u}_{i}), \qquad (A43)$$

$$u_{i+\frac{1}{2}}^{(2)} = \frac{1}{12}(\bar{u}_{i-2} - 5\bar{u}_{i-1} + 13\bar{u}_{i} + 3\bar{u}_{i+1}), \qquad (A44)$$

$$u_{i+\frac{1}{2}}^{(3)} = \frac{1}{12}(-\bar{u}_{i-1} + 7\bar{u}_{i} + 7\bar{u}_{i+1} - \bar{u}_{i+2}), \qquad (A45)$$

$$u_{i+\frac{1}{2}}^{(4)} = \frac{1}{12}(3\bar{u}_{i} + 13\bar{u}_{i+1} - 5\bar{u}_{i+2} + \bar{u}_{i+3}). \qquad (A46)$$

In the same way as before, the corresponding interpolated values for the point $x_{i-\frac{1}{2\sqrt{3}}}$ are

$$u_{i-\frac{1}{2\sqrt{3}}}^{(1)} = \frac{1}{216}(+11\sqrt{3}\bar{u}_{i-3} - 51\sqrt{3}\bar{u}_{i-2} + 105\sqrt{3}\bar{u}_{i-1} + 216\bar{u}_{i} - 65\sqrt{3}\bar{u}_{i}), \qquad (A47)$$

$$u_{i-\frac{1}{2\sqrt{3}}}^{(2)} = \frac{1}{216}(-7\sqrt{3}\bar{u}_{i-2} + 39\sqrt{3}\bar{u}_{i-1} + 216\bar{u}_{i} - 21\sqrt{3}\bar{u}_{i} - 11\sqrt{3}\bar{u}_{i+1}), \qquad (A48)$$

$$u_{i-\frac{1}{2\sqrt{3}}}^{(3)} = \frac{1}{216}(+11\sqrt{3}\bar{u}_{i-1} + 216\bar{u}_{i} + 21\sqrt{3}\bar{u}_{i} - 39\sqrt{3}\bar{u}_{i+1} + 7\sqrt{3}\bar{u}_{i+2}), \qquad (A49)$$

$$u_{i-\frac{1}{2\sqrt{3}}}^{(4)} = \frac{1}{216}(+216\bar{u}_{i} + 65\sqrt{3}\bar{u}_{i} - 105\sqrt{3}\bar{u}_{i+1} + 51\sqrt{3}\bar{u}_{i+2} - 11\sqrt{3}\bar{u}_{i+3}), \qquad (A50)$$
and the corresponding interpolated values for the point $x_{i+\frac{1}{2\sqrt{3}}}^{}$ are

$$ u_{i+\frac{1}{2\sqrt{3}}}^{(1)} = \frac{1}{216} \left( -11\sqrt{3}\bar{u}_{i-3} + 51\sqrt{3}\bar{u}_{i-2} - 105\sqrt{3}\bar{u}_{i-1} + 216\bar{u}_i + 65\sqrt{3}\bar{u}_i \right), \quad (A51) $$

$$ u_{i+\frac{1}{2\sqrt{3}}}^{(2)} = \frac{1}{216} \left( +7\sqrt{3}\bar{u}_{i-2} - 39\sqrt{3}\bar{u}_{i-1} + 216\bar{u}_i + 21\sqrt{3}\bar{u}_i + 11\sqrt{3}\bar{u}_{i+1} \right), \quad (A52) $$

$$ u_{i+\frac{1}{2\sqrt{3}}}^{(3)} = \frac{1}{216} \left( -11\sqrt{3}\bar{u}_{i-1} + 216\bar{u}_i - 21\sqrt{3}\bar{u}_i + 39\sqrt{3}\bar{u}_{i+1} - 7\sqrt{3}\bar{u}_{i+2} \right), \quad (A53) $$

$$ u_{i+\frac{1}{2\sqrt{3}}}^{(4)} = \frac{1}{216} \left( +216\bar{u}_i - 65\sqrt{3}\bar{u}_i + 105\sqrt{3}\bar{u}_{i+1} - 51\sqrt{3}\bar{u}_{i+2} + 11\sqrt{3}\bar{u}_{i+3} \right). \quad (A54) $$

If we choose the large stencil

$$ S = \{I_{i-3}, I_{i-2}, I_{i-1}, I_i, I_{i+1}, I_{i+2}, I_{i+3}\} \quad (A55) $$

which is the union of all 4 stencils $S_k$, then we are able to find an interpolation polynomial $p(x)$ of degree at most 6, satisfying $p(x_j) = u_j$ for $i-N \le j \le i+N$ and giving the approximations $u_{i-\frac{1}{2}} \equiv p(x_{i-\frac{1}{2}})$ and $u_{i+\frac{1}{2}} \equiv p(x_{i+\frac{1}{2}})$

$$ u_{i-\frac{1}{2}} = \frac{1}{420} \left( +4\bar{u}_{i-3} - 38\bar{u}_{i-2} + 214\bar{u}_{i-1} + 319\bar{u}_i - 101\bar{u}_{i+1} + 25\bar{u}_{i+2} - 3\bar{u}_{i+3} \right), \quad (A56) $$

$$ u_{i+\frac{1}{2}} = \frac{1}{420} \left( -3\bar{u}_{i-3} + 25\bar{u}_{i-2} - 101\bar{u}_{i-1} + 319\bar{u}_i + 214\bar{u}_{i+1} - 38\bar{u}_{i+2} + 4\bar{u}_{i+3} \right), \quad (A57) $$

and the approximations $u_{i-\frac{1}{2\sqrt{3}}} \equiv p(x_{i-\frac{1}{2\sqrt{3}}})$ and $u_{i+\frac{1}{2\sqrt{3}}} \equiv p(x_{i+\frac{1}{2\sqrt{3}}})$

$$ u_{i-\frac{1}{2\sqrt{3}}} = \left( +3717\sqrt{3}\bar{u}_{i-3} + 50\bar{u}_{i-3} - 32508\sqrt{3}\bar{u}_{i-2} - 552\bar{u}_{i-2} + 144585\sqrt{3}\bar{u}_{i-1} + 1758\bar{u}_{i-1} + 1086128\bar{u}_i - 144585\sqrt{3}\bar{u}_{i+1} + 1758\bar{u}_{i+1} + 32508\sqrt{3}\bar{u}_{i+2} - 552\bar{u}_{i+2} - 3717\sqrt{3}\bar{u}_{i+3} + 50\bar{u}_{i+3} \right) / 1088640, \quad (A58) $$

$$ u_{i+\frac{1}{2\sqrt{3}}} = \left( -3717\sqrt{3}\bar{u}_{i-3} + 50\bar{u}_{i-3} + 32508\sqrt{3}\bar{u}_{i-2} - 552\bar{u}_{i-2} - 144585\sqrt{3}\bar{u}_{i-1} + 1758\bar{u}_{i-1} + 1086128\bar{u}_i + 144585\sqrt{3}\bar{u}_{i+1} + 1758\bar{u}_{i+1} - 32508\sqrt{3}\bar{u}_{i+2} - 552\bar{u}_{i+2} + 3717\sqrt{3}\bar{u}_{i+3} + 50\bar{u}_{i+3} \right) / 1088640, \quad (A59) $$

provided that the function is smooth in the large stencil S.

## A5 WENO algorithm

The WENO idea is to choose the final approximation as a convex combination of the approximations $u_{i+\frac{1}{2}}^{(k)}$

$$ u_{i+\frac{1}{2}} = \sum_{k=1}^{N+1} w_k u_{i+\frac{1}{2}}^{(k)}, \quad (A60) $$

where $w_k \ge 0$, $\sum_{k=1}^{N+1} w_k = 1$.

<table><caption>Table A1. Linear weights for the WENO3, WENO5, and WENO7 schemes.</caption><thead><tr><th>N</th><th>γ</th><th>x<sub>i-1/2</sub></th><th>x<sub>i+1/2</sub></th><th>x<sub>i-1/2√3</sub></th><th>x<sub>i+1/2</sub></th></tr></thead><tbody><tr><td rowspan="2">1</td><td>γ<sub>1</sub></td><td><sup>2</sup>&frasl;<sub>3</sub></td><td><sup>1</sup>&frasl;<sub>3</sub></td><td><sup>1</sup>&frasl;<sub>2</sub></td><td><sup>1</sup>&frasl;<sub>2</sub></td></tr><tr><td>γ<sub>2</sub></td><td><sup>1</sup>&frasl;<sub>3</sub></td><td><sup>2</sup>&frasl;<sub>3</sub></td><td><sup>1</sup>&frasl;<sub>2</sub></td><td><sup>1</sup>&frasl;<sub>2</sub></td></tr><tr><td rowspan="3">2</td><td>γ<sub>1</sub></td><td><sup>3</sup>&frasl;<sub>10</sub></td><td><sup>1</sup>&frasl;<sub>10</sub></td><td><sup>210+√3</sup>&frasl;<sub>1080</sub></td><td><sup>210-√3</sup>&frasl;<sub>1080</sub></td></tr><tr><td>γ<sub>2</sub></td><td><sup>3</sup>&frasl;<sub>5</sub></td><td><sup>3</sup>&frasl;<sub>5</sub></td><td><sup>11</sup>&frasl;<sub>18</sub></td><td><sup>11</sup>&frasl;<sub>18</sub></td></tr><tr><td>γ<sub>3</sub></td><td><sup>1</sup>&frasl;<sub>10</sub></td><td><sup>3</sup>&frasl;<sub>10</sub></td><td><sup>210-√3</sup>&frasl;<sub>1080</sub></td><td><sup>210+√3</sup>&frasl;<sub>1080</sub></td></tr><tr><td rowspan="4">3</td><td>γ<sub>1</sub></td><td><sup>4</sup>&frasl;<sub>35</sub></td><td><sup>1</sup>&frasl;<sub>35</sub></td><td><sup>11151+50√3</sup>&frasl;<sub>166320</sub></td><td><sup>11151-50√3</sup>&frasl;<sub>166320</sub></td></tr><tr><td>γ<sub>2</sub></td><td><sup>18</sup>&frasl;<sub>35</sub></td><td><sup>12</sup>&frasl;<sub>35</sub></td><td><sup>168021+1174√3</sup>&frasl;<sub>388080</sub></td><td><sup>168021-1174√3</sup>&frasl;<sub>388080</sub></td></tr><tr><td>γ<sub>3</sub></td><td><sup>12</sup>&frasl;<sub>35</sub></td><td><sup>18</sup>&frasl;<sub>35</sub></td><td><sup>168021-1174√3</sup>&frasl;<sub>388080</sub></td><td><sup>168021+1174√3</sup>&frasl;<sub>388080</sub></td></tr><tr><td>γ<sub>4</sub></td><td><sup>1</sup>&frasl;<sub>35</sub></td><td><sup>4</sup>&frasl;<sub>35</sub></td><td><sup>11151-50√3</sup>&frasl;<sub>166320</sub></td><td><sup>11151+50√3</sup>&frasl;<sub>166320</sub></td></tr></tbody></table>

## A6 Linear weights

The approximations calculated by using the large stencils can be written as a linear convex combination of the approximations $u_{i+\frac{1}{2}}^{(k)}$ based on the small stencils $S_i$

$$ u_{i+\frac{1}{2}} = \sum_{k=1}^{N+1} \gamma_k u_{i+\frac{1}{2}}^{(k)} \quad (A61) $$

where the constants $\gamma_j$ satisfy $\sum_{j=1}^{N+1} \gamma_j = 1$, are usually referred to as the linear weights in the WENO literature. We provide the linear weights for the points of interest considered before, namely, $x_{i-\frac{1}{2}}$, $x_{i+\frac{1}{2}}$, $x_{i-\frac{1}{2\sqrt{3}}}$, and $x_{i+\frac{1}{2\sqrt{3}}}$, in the table A1.

To find the linear weights, simply solve for $\gamma_j$ in the following equation system

$$ \sum_{i=1}^{N+1} \gamma_i p^{(i)}(x) = \sum_{k=i-N}^{i+N} \alpha_k u_k(x), \quad (A62) $$

where the right hand side was determined with interpolation in the large stencil. This system is not always solvable, and for some cases, the linear weights are negative (see Shi et al. (2002)). All this means that the WENO reconstruction procedure can only be applied for certain points in the interval $I_i = (x_{i-\frac{1}{2}}, x_{i+\frac{1}{2}})$.

## A7 Nonlinear weights

The nonlinear weights $w_k$ satisfy the following requirements

* $\omega_k \approx \gamma_k$ if $u(x)$ is smooth in the big stencil S;

* $\omega_k \approx 0$ if $u(x)$ has a discontinuity in the stencil $S_k$ but is small in at least one of the other stencils.

It can be shown (Jiang & Shu 1996) that, as long as

$$ \omega_k = \gamma_k + O(\Delta x^N), \quad (A63) $$

the WENO interpolation $u_{i+\frac{1}{2}}$ is $(2N+1)$th order accurate,

$$ u_{i+\frac{1}{2}} - u(x_{i+\frac{1}{2}}) = O(\Delta x^{2N+1}), \quad (A64) $$

when the function $u(x)$ is smooth in the large stencil S. The second requirement above would guarantee a nonoscillatory, at least Nth order accurate WENO approximation $u_{i+\frac{1}{2}}$, since the contribution from any stencil containing the discontinuity of $u(x)$ has an essentially zero weight.
The choice of the nonlinear weights $\omega_k$ relies on the smoothness indicator $\beta_k$, which measures the relative smoothness of the function $u(x)$ in the stencil $S_k$. The larger this smoothness indicator $\beta_k$, the less smooth the function $u(x)$ is in the stencil $S_k$. In the most of the WENO papers, the smoothness indicator is chosen as in (Jiang & Shu 1996)

$$ \beta_k = \sum_{l=1}^{N} \Delta x^{2l-1} \int_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}} \left( \frac{d^l}{dx^l} p_k(x) \right)^2 dx, \quad (A65) $$

where $N$ is the polynomial degree of $p_k(x)$ (in our examples, $N = 1, 2, 3$). This is clearly just a scaled sum of the square $L^2$ norms of all the derivatives of the relevant interpolation polynomial $p_k(x)$ in the relevant interval $[x_{i-\frac{1}{2}}, x_{i+\frac{1}{2}}]$, where the interpolating point is located. The scaling factor $\Delta x^{2l-1}$ is to make sure that the final explicit formulae for the smoothness indicators do not depend on the mesh size $\Delta x$. Let us write the smoothness indicators for polynomials of degree $N = 1$

$$ \beta_1 = (\bar{u}_{i-1} - \bar{u}_i)^2, \quad (A66) $$

$$ \beta_2 = (\bar{u}_i - \bar{u}_{i+1})^2, \quad (A67) $$

and for polynomials of degree $N = 2$ we get

$$ \beta_1 = \frac{1}{3} \left( \begin{aligned} & + 4\bar{u}_{i-2}^2 - 19\bar{u}_{i-1}\bar{u}_{i-2} + 11\bar{u}_i\bar{u}_{i-2} \\ & + 25\bar{u}_{i-1}^2 + 10\bar{u}_i^2 - 31\bar{u}_{i-1}\bar{u}_i \end{aligned} \right), \quad (A68) $$

$$ \beta_2 = \frac{1}{3} \left( \begin{aligned} & + 4\bar{u}_{i-1}^2 - 13\bar{u}_i\bar{u}_{i-1} + 5\bar{u}_{i+1}\bar{u}_{i-1} \\ & + 13\bar{u}_i^2 + 4\bar{u}_{i+1}^2 - 13\bar{u}_i\bar{u}_{i+1} \end{aligned} \right), \quad (A69) $$

$$ \beta_3 = \frac{1}{3} \left( \begin{aligned} & + 10\bar{u}_i^2 - 31\bar{u}_{i+1}\bar{u}_i + 11\bar{u}_{i+2}\bar{u}_i \\ & + 25\bar{u}_{i+1}^2 + 4\bar{u}_{i+2}^2 - 19\bar{u}_{i+1}\bar{u}_{i+2} \end{aligned} \right), \quad (A70) $$

and for polynomials of degree $N = 3$ we get

$$ \beta_1 = \frac{1}{240} \left( \begin{aligned} & + 547\bar{u}_{i-3}^2 - 3882\bar{u}_{i-3}\bar{u}_{i-2} + 4642\bar{u}_{i-3}\bar{u}_{i-1} \\ & - 1854\bar{u}_{i-3}\bar{u}_i + 7043\bar{u}_{i-2}^2 - 17246\bar{u}_{i-2}\bar{u}_{i-1} \\ & + 7042\bar{u}_{i-2}\bar{u}_i + 11003\bar{u}_{i-1}^2 - 9402\bar{u}_{i-1}\bar{u}_i \\ & + 2107\bar{u}_i^2 \end{aligned} \right), \quad (A71) $$

$$ \beta_2 = \frac{1}{240} \left( \begin{aligned} & + 267\bar{u}_{i-2}^2 - 1642\bar{u}_{i-2}\bar{u}_{i-1} + 1602\bar{u}_{i-2}\bar{u}_i \\ & - 494\bar{u}_{i-2}\bar{u}_{i+1} + 2843\bar{u}_{i-1}^2 - 5966\bar{u}_{i-1}\bar{u}_i \\ & + 1922\bar{u}_{i-1}\bar{u}_{i+1} + 3443\bar{u}_i^2 - 2522\bar{u}_i\bar{u}_{i+1} \\ & + 547\bar{u}_{i+1}^2 \end{aligned} \right), \quad (A72) $$

$$ \beta_3 = \frac{1}{240} \left( \begin{aligned} & + 547\bar{u}_{i-1}^2 - 2522\bar{u}_{i-1}\bar{u}_i + 1922\bar{u}_{i-1}\bar{u}_{i+1} \\ & - 494\bar{u}_{i-1}\bar{u}_{i+2} + 3443\bar{u}_i^2 - 5966\bar{u}_i\bar{u}_{i+1} \\ & + 1602\bar{u}_i\bar{u}_{i+2} + 2843\bar{u}_{i+1}^2 - 1642\bar{u}_{i+1}\bar{u}_{i+2} \\ & + 267\bar{u}_{i+2}^2 \end{aligned} \right), \quad (A73) $$

$$ \beta_4 = \frac{1}{240} \left( \begin{aligned} & + 2107\bar{u}_i^2 - 9402\bar{u}_i\bar{u}_{i+1} + 7042\bar{u}_i\bar{u}_{i+2} \\ & - 1854\bar{u}_i\bar{u}_{i+3} + 11003\bar{u}_{i+1}^2 - 17246\bar{u}_{i+1}\bar{u}_{i+2} \\ & + 4642\bar{u}_{i+1}\bar{u}_{i+3} + 7043\bar{u}_{i+2}^2 - 3882\bar{u}_{i+2}\bar{u}_{i+3} \\ & + 547\bar{u}_{i+3}^2 \end{aligned} \right), \quad (A74) $$

Notice that these smoothness indicators are quadratic functions of the values of $u(x)$ in the relevant stencils. Equipped with these

smoothness indicators, we can now define the nonlinear weights as

$$ \omega_k = \frac{\tilde{\omega}_k}{\sum_{i=1}^{N} \tilde{\omega}_i}, \quad \text{with} \quad \tilde{\omega}_k = \frac{\gamma_k}{(\varepsilon + \beta_k)^2}. \quad (A75) $$

Here $\varepsilon$ is a small positive number used to avoid the denominator becoming zero and is typically chosen to be $\varepsilon = 10^{-6}$, but in our calculations we have used $\varepsilon = 10^{-24}$.

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