# A $4^{\rm th}$-order accurate finite volume method for ideal classical and
special relativistic MHD based on pointwise reconstructions

V. Berta  A. Mignone  M. Bugli  G. Mattia  Dipartimento di Fisica, Università
di Torino, via Pietro Giuria 1, 10125 Torino, Italy  Université Paris-Saclay,
Université Paris Cité, CEA, CNRS, AIM, 91191, Gif-sur-Yvette, France  INFN,
Sezione di Firenze, Via G. Sansone 1, I-50019 Sesto Fiorentino, Italy

###### Abstract

We present a novel implementation of a genuinely $4^{\rm th}$-order accurate
finite volume scheme for multidimensional classical and special relativistic
magnetohydrodynamics (MHD) based on the constrained transport (CT) formalism.
The scheme introduces several novel aspects when compared to its predecessors
yielding a more efficient computational tool. Among the most relevant ones,
our scheme exploits pointwise to pointwise reconstructions (rather than one-
dimensional finite volume ones), employs the generic upwind constrained
transport averaging and sophisticated limiting strategies that include both a
discontinuity detector and an order reduction procedure. Selected numerical
benchmarks demonstrate the accuracy and robustness of the method.

###### keywords:

Magnetohydrodynamics (MHD) , high-order finite volume methods , constrained
transport , Riemann solvers

††journal: Journal of Computational Physics

##  1 Introduction

The modeling of astrophysical plasmas is nowadays demanding for more efficient
and accurate numerical schemes for solving the equations of classical and
special relativistic magnetohydrodynamics (MHD, RMHD). Physical interest is
increasingly directed towards the investigation of highly nonlinear flows
featuring both smooth and discontinuous solutions, with a great deal of
attention being devoted to enhance the accuracy of the underlying numerical
methods while retaining stability and computational efficiency.

Traditionally, in the context of finite volume (FV) or finite difference (FD)
schemes, a great variety of $2^{\rm nd}$-order methods have been developed to
solve the MHD and RMHD equations (e.g., [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11,
12, 13, 14, 15, 16, 17, 18]). However, in the last decade several efforts have
been made to design higher than $2^{\rm nd}$-order numerical methods both for
FD (e.g., [19, 20, 15, 21, 22, 23, 24]) as well as for FV (e.g., [25, 26, 27,
28, 29, 30, 31, 32]) approaches, leading to the dawn of high-order numerical
codes (e.g. [33, 34, 35, 36, 37, 38]) in plasma physics and computational
fluid-dynamics. Most $4^{\rm th}$-order schemes, in fact, can reach
unprecedented performances overstepping the limits of traditional $2^{\rm
nd}$-order frameworks because of their intrinsic lower dissipation properties.
In addition, smooth solutions are improved at much faster rates, yielding
enhanced convergence and henceforth a computational efficiency that increases
with the dimensionality of the problem and with resolution.

In this paper we present a genuine $4^{\rm th}$-order accurate finite volume
method for the numerical solution of the ideal classical and special
relativistic MHD equations. Such method draws on the work originally
introduced by McCorquodale & Colella [28] and then refined by Felker & Stone
[36], but with a number of important improvements and differences. In fact,
our numerical scheme is introducing, for the first time in the literature of
FV schemes, a version of the MP5 [39] and the WENOZ [40, 41, 42] spatial
reconstruction algorithms based on pointwise values, which are instead
traditionally used in FD schemes. As we shall demonstrate, pointwise
reconstructions lower the algorithmic complexity of the overall scheme by
retaining the $4^{\rm th}$-order accuracy with fewer operations per step and
without even enlarging the stencil size. Moreover, this new feature reduces
the number of Riemann problems to be solved per direction, leading to a more
efficient and cost-effective scheme. Concurrently, robustness is ensured by
means of a local discontinuity detector to distinguish smooth solutions from
discontinuous ones when determining the function local point value from its
cell average. Regarding temporal integration, we employ a semi-discrete
approach based on a five-stage $4^{\rm th}$-order explicit strong stability
preserving Runge-Kutta method (eSSPRK(5,4), [43, 44]). The solenoidal
condition of the magnetic field is fulfilled by means of a high-order
formulation of the generalized upwind constrained transport (UCT, [45, 3, 5,
46, 47]) algorithm, which ensures a divergence-free magnetic field up to
machine accuracy, and works well with any generic upwind average [18].

The paper is structured as follow: in §2 we will briefly describe the MHD and
RMHD equations and how they are discretized in a general FV framework together
with the notation used throughout this work. In §3 an extensive description of
the novel numerical method will be provided. In §4 we will illustrate our
method to control the solenoidal condition. Numerical tests and conclusions
will be given, respectively, in §5 and §6.

##  2 The ideal MHD and RMHD equations

The ideal MHD equations can be split out in two coupled sub-systems [19, 47,
18], the first one being a time-dependent hyperbolic system for the evolution
of the conservative flow variables (the mass density $\rho$, the momentum
density $\rho\mathbf{v}$, and the energy density $\altmathcal{E}$). This sub-
system can be written in compact form by introducing a state vector of
conservative quantities, $U$, and a rank-2 tensor, $\mathsf{F}$, whose rows
represent the fluxes relative to each component of $U$

$\frac{\partial U}{\partial t}+\nabla\cdot\mathsf{F}=0\,,$ (1)

where

$U=\left(\begin{array}[]{l}\rho\\\ \vskip 6.0pt plus 2.0pt minus
2.0pt\cr\rho\mathbf{v}\\\ \vskip 6.0pt plus 2.0pt minus
2.0pt\cr\altmathcal{E}\end{array}\right)\,,\quad\mathsf{F}=\left(\begin{array}[]{c}\rho\mathbf{v}\\\
\vskip 6.0pt plus 2.0pt minus
2.0pt\cr\rho\mathbf{v}\mathbf{v}-\mathbf{B}\mathbf{B}+\mathsf{I}p_{t}\\\
\vskip 6.0pt plus 2.0pt minus
2.0pt\cr(\altmathcal{E}+p_{t})\mathbf{v}-(\mathbf{v}\cdot\mathbf{B})\mathbf{B}\end{array}\right)^{\intercal}\,.\quad\vspace{0.3
cm}$ (2)

In Eq. (2), $\mathbf{v}=(v_{x},\,v_{y},\,v_{z})$ is the fluid velocity,
$\mathbf{B}=(B_{x},\,B_{y},\,B_{z})$ is the magnetic field, $p_{t}=p+B^{2}/2$
is the total pressure expressed as the sum of the thermal and magnetic
pressure, and $\mathsf{I}$ is the identity matrix. The total energy density is
composed by a kinetic, a thermal and a magnetic term

$\altmathcal{E}=\frac{1}{2}\rho\mathbf{v}^{2}+\frac{p}{\Gamma-1}+\frac{\mathbf{B}^{2}}{2}\,,$
(3)

while $\Gamma$ is the specific heat ratio for an adiabatic equation of state.
On the other hand, the second sub-system consists of the induction equation
for the evolution of the magnetic field

$\frac{\partial\mathbf{B}}{\partial t}+c\nabla\times\mathbf{E}=0.$ (4)

Due to the presence of the curl operator instead of the divergence, the latter
equation cannot be considered a simple extension of the HD subsystem [19].
Additionally, in the ideal MHD framework, the electric field $\mathbf{E}$ is a
function of the fluid velocity $\mathbf{v}$ and the magnetic field
$\mathbf{B}$

$c\mathbf{E}=-\mathbf{v}\times\mathbf{B}\,,$ (5)

where $c$ is the speed of light.

Analytically, the commutativity of spatial derivatives endows the induction
equation with the solenoidal condition for the magnetic field

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

This stationary condition, once satisﬁed at $t=0$, is analytically preserved
at $t>0$. On the other hand, since numerical derivatives do not commute, the
condition expressed in Eq. (6) is not automatically fulfilled by numerical
schemes, leading to the rising of magnetic monopoles which cause unphysical
plasma transport and dissipation of momentum and energy in ideal frameworks
[48]. For these reasons, HD numerical schemes cannot be straightforwardly
extended to MHD systems without introducing some corrections to numerically
constraint Eq. (6) to be satisfied.

Similarly to the MHD case, the ideal RMHD equations can be expressed in the
form given by Eq. (1) with

$U=\left(\begin{array}[]{l}D\\\ \vskip 6.0pt plus 2.0pt minus
2.0pt\cr\mathbf{m}\\\ \vskip 6.0pt plus 2.0pt minus
2.0pt\cr\altmathcal{E}\end{array}\right)\,,\quad\mathsf{F}=\left(\begin{array}[]{c}D\mathbf{v}\\\
\vskip 6.0pt plus 2.0pt minus 2.0pt\cr
w_{t}\gamma^{2}\mathbf{v}\mathbf{v}-\mathbf{b}\mathbf{b}+\mathsf{I}p_{t}\\\
\vskip 6.0pt plus 2.0pt minus
2.0pt\cr\mathbf{m}\end{array}\right)^{\intercal}\,,\quad\vspace{0.3 cm}$ (7)

where $D=\gamma\rho$, $\mathbf{m}=(\rho
h\gamma^{2}+|\mathbf{B}|^{2})\mathbf{v}-(\mathbf{v}\cdot\mathbf{B})\mathbf{B}$,
and $\altmathcal{E}=\rho
h\gamma^{2}-p+|\mathbf{B}|^{2}/2+(|\mathbf{v}|^{2}|\mathbf{B}|^{2}-(\mathbf{v}\cdot\mathbf{B})^{2})/2$
are respectively the relativistic mass density, momentum density, and energy
density, while $w_{t}=\rho
h+|\mathbf{B}|^{2}/\gamma^{2}+(\mathbf{v}\cdot\mathbf{B})^{2}$ and
$\mathbf{b}=\mathbf{B}/\gamma+\gamma(\mathbf{v}\cdot\mathbf{B})\mathbf{v}$. In
the previous expressions we introduced the specific enthalpy $h$, the Lorentz
factor $\gamma=(1-|\mathbf{v}|^{2})^{-1/2}$ and we set the speed of light $c$
to unity. In addition, the same induction equation as in the MHD regime needs
to be added to Eq. (7), thus requiring the same numerical scheme to control
the divergence-free condition shown in Eq. (6).

###  2.1 Notations

We employ a Cartesian coordinate system with unit vectors
$\hat{\mathbf{e}}_{x}=(1,0,0)$, $\hat{\mathbf{e}}_{y}=(0,1,0)$, and
$\hat{\mathbf{e}}_{z}=(0,0,1)$ and mantain the same notation formalism adopted
by Mignone & Del Zanna [18]. In each direction the coordinate spacing $\Delta
x$, $\Delta y$, and $\Delta z$ is uniform. Computational cells are centered at
$(x_{i},\,y_{j},\,z_{k})$ and delimited by the six interfaces orthogonal to
the coordinate axis centered, respectively, at
$(x_{i\pm\frac{1}{2}},\,y_{j},\,z_{k})$,
$(x_{i},\,y_{j\pm\frac{1}{2}},\,z_{k})$, and
$(x_{i},\,y_{j},\,z_{k\pm\frac{1}{2}})$.

According to Gauss’ theorem, the conserved quantities are evolved as volume
averages $\left<{U}\right>_{{\boldsymbol{c}}}$ over the cell volume

$\left<{U}\right>_{{\boldsymbol{c}}}\equiv\frac{1}{\Delta x\Delta y\Delta
z}\int U(x,y,z,t)\,dxdydz\,,$ (8)

where the ${\boldsymbol{c}}$ subscript is a shorthand notation for $(i,j,k)$,
by means of the neat difference of the surface-averaged representation of the
fluxes at zone interfaces:

$\begin{array}[]{l}\displaystyle\hat{F}_{x,{\mathbf{x}_{f}}}\equiv\frac{1}{\Delta
y\Delta
z}\int\hat{\mathbf{e}}_{x}\cdot\mathsf{F}\big{(}U(x_{i+\frac{1}{2}},y,z,t)\big{)}\,dy\,dz\,,\\\
\vskip 6.0pt plus 2.0pt minus
2.0pt\cr\displaystyle\hat{F}_{y,{\mathbf{y}_{f}}}\equiv\frac{1}{\Delta x\Delta
z}\int\hat{\mathbf{e}}_{y}\cdot\mathsf{F}\big{(}U(x,y_{j+\frac{1}{2}},z,t)\big{)}\,dz\,dx\,,\\\
\vskip 6.0pt plus 2.0pt minus
2.0pt\cr\displaystyle\hat{F}_{z,{\mathbf{z}_{f}}}\equiv\frac{1}{\Delta x\Delta
y}\int\hat{\mathbf{e}}_{z}\cdot\mathsf{F}\big{(}U(x,y,z_{k+\frac{1}{2}},t)\big{)}\,dx\,dy\,.\end{array}$
(9)

In the finite-volume fashion, a semi-discrete method of lines approximates the
partial differential equation (PDE) associated with Eq. (1), yielding an
ordinary differential equation (ODE) in the time variable

$\displaystyle\frac{d\left<{U}\right>_{{\boldsymbol{c}}}}{dt}=-\left(\frac{\hat{F}_{x,{\mathbf{x}_{f}}}-\hat{F}_{x,{\mathbf{x}_{f}}-\hat{\mathbf{e}}_{x}}}{\Delta
x}+\frac{\hat{F}_{y,{\mathbf{y}_{f}}}-\hat{F}_{y,{\mathbf{y}_{f}}-\hat{\mathbf{e}}_{y}}}{\Delta
y}+\frac{\hat{F}_{z,{\mathbf{z}_{f}}}-\hat{F}_{z,{\mathbf{z}_{f}}-\hat{\mathbf{e}}_{z}}}{\Delta
z}\right)\,.$ (10)

Eqns. (8-10) frame the discretization for zone-centered variables. In
constrained transport MHD instead, the magnetic field update relies on a
discrete version of Stoke’s theorem [49, 50, 45, 5], where the magnetic field
is treated as a staggered primary variable. In this representation, each
component of the magnetic field

$\mathbf{\hat{B}}_{f}\equiv\left(\begin{array}[]{l}\hat{B}_{x,i+\frac{1}{2},j,k}\\\
\vskip 6.0pt plus 2.0pt minus 2.0pt\cr\hat{B}_{y,i,j+\frac{1}{2},k}\\\ \vskip
6.0pt plus 2.0pt minus
2.0pt\cr\hat{B}_{z,i,j,k+\frac{1}{2}}\end{array}\right)=\left(\begin{array}[]{l}\hat{B}_{x,{\mathbf{x}_{f}}}\\\
\vskip 6.0pt plus 2.0pt minus 2.0pt\cr\hat{B}_{y,{\mathbf{y}_{f}}}\\\ \vskip
6.0pt plus 2.0pt minus
2.0pt\cr\hat{B}_{z,{\mathbf{z}_{f}}}\end{array}\right)\,$ (11)

is evolved as an area-weighted average on the face orthogonal to the field
component. The subscripts ${\mathbf{x}_{f}}$, ${\mathbf{y}_{f}}$ and
${\mathbf{z}_{f}}$ identify the location of the faces, i.e.,
${\mathbf{x}_{f}}\equiv(i+\frac{1}{2},j,k)$,
${\mathbf{y}_{f}}\equiv(i,j+\frac{1}{2},k)$, and
${\mathbf{z}_{f}}\equiv(i,j,k+\frac{1}{2})$.

Furthermore, Eq. (4) prescribes that the variation in time of the magnetic
flux is equal to the line integral of the electric field
$\bar{\mathbf{E}}_{e}$ over zone edges, namely, the electromotive force (EMF)

$\bar{\mathbf{E}}_{e}\equiv\left(\begin{array}[]{l}\bar{E}_{x,i,j+\frac{1}{2},k+\frac{1}{2}}\\\
\vskip 6.0pt plus 2.0pt minus
2.0pt\cr\bar{E}_{y,i+\frac{1}{2},j,k+\frac{1}{2}}\\\ \vskip 6.0pt plus 2.0pt
minus
2.0pt\cr\bar{E}_{z,i+\frac{1}{2},j+\frac{1}{2},k}\end{array}\right)=\left(\begin{array}[]{l}\bar{E}_{x,{\mathbf{x}_{e}}}\\\
\vskip 6.0pt plus 2.0pt minus 2.0pt\cr\bar{E}_{y,{\mathbf{y}_{e}}}\\\ \vskip
6.0pt plus 2.0pt minus
2.0pt\cr\bar{E}_{z,{\mathbf{z}_{e}}}\end{array}\right)\,,$ (12)

where:

$\begin{array}[]{l}\displaystyle\bar{E}_{x,{\mathbf{x}_{e}}}\equiv\frac{1}{\Delta
x}\int E_{x}(x,y_{j+\frac{1}{2}},z_{k+\frac{1}{2}},t)\,dx\,,\\\ \vskip 6.0pt
plus 2.0pt minus
2.0pt\cr\displaystyle\bar{E}_{y,{\mathbf{y}_{e}}}\equiv\frac{1}{\Delta y}\int
E_{y}(x_{i+\frac{1}{2}},y,z_{k+\frac{1}{2}},t)\,dy\,,\\\ \vskip 6.0pt plus
2.0pt minus
2.0pt\cr\displaystyle\bar{E}_{z,{\mathbf{z}_{e}}}\equiv\frac{1}{\Delta z}\int
E_{z}(x_{i+\frac{1}{2}},y_{j+\frac{1}{2}},z,t)\,dz\,.\end{array}$ (13)

Similarly to the magnetic field, the edge-centered electric field positions
are labeled as ${\mathbf{x}_{e}}\equiv(i,j+\frac{1}{2},k+\frac{1}{2})$,
${\mathbf{y}_{e}}\equiv(i+\frac{1}{2},j,k+\frac{1}{2})$, and
${\mathbf{z}_{e}}\equiv(i+\frac{1}{2},j+\frac{1}{2},k)$. Given all these
elements, it is possible to retrieve the discretized version of Stokes theorem

$\begin{array}[]{lcl}\displaystyle\frac{d\hat{B}_{x,{\mathbf{x}_{f}}}}{dt}&=&\displaystyle-\left(\frac{\bar{E}_{z,{\mathbf{z}_{e}}}-\bar{E}_{z,{\mathbf{z}_{e}}-\hat{\mathbf{e}}_{y}}}{\Delta
y}-\frac{\bar{E}_{y,{\mathbf{y}_{e}}}-\bar{E}_{y,{\mathbf{y}_{e}}-\hat{\mathbf{e}}_{z}}}{\Delta
z}\right)\,,\\\ \vskip 6.0pt plus 2.0pt minus
2.0pt\cr\displaystyle\frac{d\hat{B}_{y,{\mathbf{y}_{f}}}}{dt}&=&\displaystyle-\left(\frac{\bar{E}_{x,{\mathbf{x}_{e}}}-\bar{E}_{x,{\mathbf{x}_{e}}-\hat{\mathbf{e}}_{z}}}{\Delta
z}-\frac{\bar{E}_{z,{\mathbf{z}_{e}}}-\bar{E}_{z,{\mathbf{z}_{e}}-\hat{\mathbf{e}}_{x}}}{\Delta
x}\right)\,,\\\ \vskip 6.0pt plus 2.0pt minus
2.0pt\cr\displaystyle\frac{d\hat{B}_{z,{\mathbf{z}_{f}}}}{dt}&=&\displaystyle-\left(\frac{\bar{E}_{y,{\mathbf{y}_{e}}}-\bar{E}_{y,{\mathbf{y}_{e}}-\hat{\mathbf{e}}_{x}}}{\Delta
x}-\frac{\bar{E}_{x,{\mathbf{x}_{e}}}-\bar{E}_{x,{\mathbf{x}_{e}}-\hat{\mathbf{e}}_{y}}}{\Delta
y}\right)\,,\end{array}$ (14)

where we set $c=1$, thus verifying

$\frac{d}{dt}\left(\frac{\hat{B}_{x,{\mathbf{x}_{f}}}-\hat{B}_{x,{\mathbf{x}_{f}}-\hat{\mathbf{e}}_{x}}}{\Delta
x}+\frac{\hat{B}_{y,{\mathbf{y}_{f}}}-\hat{B}_{y,{\mathbf{y}_{f}}-\hat{\mathbf{e}}_{y}}}{\Delta
y}+\frac{\hat{B}_{z,{\mathbf{z}_{f}}}-\hat{B}_{z,{\mathbf{z}_{f}}-\hat{\mathbf{e}}_{z}}}{\Delta
z}\right)=0\,,$ (15)

which is valid exactly and guaranties the fulfillment of the solenoidal
condition for an initially divergence-free magnetic field. Notice that no
approximation has been made in the treatment so far. To summarize our notation
convention, given the generic physical quantity $Q$ we indicate with
$\left<{Q}\right>_{{\boldsymbol{c}}}$ its volume average. We use
$\hat{Q}_{{\mathbf{x}_{f}},{\mathbf{y}_{f}},{\mathbf{z}_{f}}}$ and
$\bar{Q}_{{\mathbf{x}_{e}},{\mathbf{y}_{e}},{\mathbf{z}_{e}}}$ to refer to the
face and edge averages, respectively. Similarly, the notations
$Q_{{\boldsymbol{c}}}$,
$Q_{{\mathbf{x}_{f}},{\mathbf{y}_{f}},{\mathbf{z}_{f}}}$, and
$Q_{{\mathbf{x}_{e}},{\mathbf{y}_{e}},{\mathbf{z}_{e}}}$ refer to the
pointwise values at the center of the cell, face, and edge, respectively.

##  3 The fourth-order accurate method

In this section we introduce the main features of our $4^{\rm th}$-order
accurate numerical scheme. Finite volume Godunov-type methods evolve the cell-
averaged conservative variables (e.g.,
$\left<{U}\right>_{\boldsymbol{c}}=\\{\left<{\rho}\right>_{\boldsymbol{c}},\left<{\rho\mathbf{v}}\right>_{\boldsymbol{c}},\left<{\altmathcal{E}}\right>_{\boldsymbol{c}}$},
while the representation of $\vec{B}$ depends on the choice of the algorithm
chosen to control the solenoidal condition and may be either a cell-centered
or a staggered quantity) and typically consist of three main stages: a
reconstruction step at cell faces, flux computation at zone interfaces via
Riemann solvers, and time integration. At $2^{\rm nd}$-order accuracy, point
values can be interchanged with volume or surface averages, so that the
transformation between conservative and primitive variables (e.g.,
$V_{c}=\\{\rho,\mathbf{v},p\\}$ same arguments as before hold for $\vec{B}$)
can be simply defined as
$V_{\boldsymbol{c}}\approx\left<{V}\right>_{{\boldsymbol{c}}}\approx{\mathcal{V}}(\left<{U}\right>_{\boldsymbol{c}})$.
This property, however, does no longer hold for higher order schemes, where
one has to distinguish the cell-centered point value from its volume average
and, likewise, the face-centered flux from its surface average. On top of
that, one-dimensional FV reconstruction schemes have to be reformulated in
order to obtain $4^{\rm th}$-order accurate left and right states at zone
interfaces. These states are best obtained using primitive or characteristic
variables in order to avoid unwanted numerical oscillations. These steps are
described in the following sections.

###  3.1 Point value recovery of primitive variables

Before reconstruction of primitive variables can take place, their values must
become available at cell centers. Since conservative schemes evolve the volume
average of conservative quantities $\left<{Q}\right>_{\boldsymbol{c}}$, their
local point value $Q_{\boldsymbol{c}}$ has to be obtained first. To $4^{\rm
th}$-order accuracy, point values can be retrieved from cell averages using
the relation

$Q_{\boldsymbol{c}}=\displaystyle\left<{Q}\right>_{\boldsymbol{c}}-\frac{\Delta\left<{Q}\right>_{\boldsymbol{c}}}{24}+O(h^{4})\,,$
(16)

and conversely

$\left<{Q}\right>_{\boldsymbol{c}}=Q_{\boldsymbol{c}}+\frac{\Delta
Q_{\boldsymbol{c}}}{24}+O(h^{4})\,,$ (17)

where Eq. (17) is effectively a high-order Simpson-like quadrature rule. In
Eqns. (16) and (17) $Q$ is a component of $U_{\boldsymbol{c}}$, while $\Delta$
is the Laplacian operator first introduced in the context of high-order
methods by McCorquodale & Colella [28], i.e.

$\Delta
Q_{{\boldsymbol{c}}}\equiv\Delta^{x}Q_{\boldsymbol{c}}+\Delta^{y}Q_{\boldsymbol{c}}+\Delta^{z}Q_{\boldsymbol{c}}\,,$
(18)

with components given by:

$\Delta^{x}Q_{{\boldsymbol{c}}}\equiv\displaystyle\left(Q_{{\boldsymbol{c}}-\hat{\mathbf{e}}_{x}}-2Q_{{\boldsymbol{c}}}+Q_{{\boldsymbol{c}}+\hat{\mathbf{e}}_{x}}\right)\,,\,\quad\Delta^{y}Q_{{\boldsymbol{c}}}\equiv\displaystyle\left(Q_{{\boldsymbol{c}}-\hat{\mathbf{e}}_{y}}-2Q_{{\boldsymbol{c}}}+Q_{{\boldsymbol{c}}+\hat{\mathbf{e}}_{y}}\right)\,,\,\quad\Delta^{z}Q_{{\boldsymbol{c}}}\equiv\displaystyle\left(Q_{{\boldsymbol{c}}-\hat{\mathbf{e}}_{z}}-2Q_{{\boldsymbol{c}}}+Q_{{\boldsymbol{c}}+\hat{\mathbf{e}}_{z}}\right)\,.$
(19)

The interested reader may refer to A for a full derivation of Eqns. (16) and
(17).

Similarly, high-order transformations need to be applied also in the
conversion between point values and face- or edge-averaged quantities, by
means of the corresponding 2D or 1D operators. At zone faces, for instance, we
use transverse Laplacian operators defined by:

$\Delta^{x}_{\bot}Q_{{\mathbf{x}_{f}}}=\Delta^{y}Q_{\mathbf{x}_{f}}+\Delta^{z}Q_{\mathbf{x}_{f}}\,,\quad\Delta^{y}_{\bot}Q_{{\mathbf{y}_{f}}}=\Delta^{x}Q_{\mathbf{y}_{f}}+\Delta^{z}Q_{\mathbf{y}_{f}}\,,\quad\Delta^{z}_{\bot}Q_{{\mathbf{z}_{f}}}=\Delta^{x}Q_{\mathbf{z}_{f}}+\Delta^{y}Q_{\mathbf{z}_{f}}\,.$
(20)

This is employed, for example, to recover the face-centered point value of the
magnetic field from its surface average:

$B_{x,{\mathbf{x}_{f}}}=\hat{B}_{x,{\mathbf{x}_{f}}}-\frac{\Delta^{x}_{\bot}\hat{B}_{x,{\mathbf{x}_{f}}}}{24},\quad
B_{y,{\mathbf{y}_{f}}}=\hat{B}_{y,{\mathbf{y}_{f}}}-\frac{\Delta^{y}_{\bot}\hat{B}_{y,{\mathbf{y}_{f}}}}{24},\quad
B_{z,{\mathbf{z}_{f}}}=\hat{B}_{z,{\mathbf{z}_{f}}}-\frac{\Delta^{z}_{\bot}\hat{B}_{z,{\mathbf{z}_{f}}}}{24}\,.$
(21)

Following Felker & Stone [36], the pointwise value of the magnetic field at
the cell center is obtained using an unlimited high-order interpolation, i.e.:

$\begin{array}[]{l}\displaystyle
B_{x,c}=-\frac{1}{16}(B_{x,{\mathbf{x}_{f}}+\hat{\mathbf{e}}_{x}}+B_{x,{\mathbf{x}_{f}}-2\hat{\mathbf{e}}_{x}})+\frac{9}{16}(B_{x,{\mathbf{x}_{f}}}+B_{x,{\mathbf{x}_{f}}-\hat{\mathbf{e}}_{x}})\,,\\\
\vskip 6.0pt plus 2.0pt minus 2.0pt\cr\displaystyle
B_{y,c}=-\frac{1}{16}(B_{y,{\mathbf{y}_{f}}+\hat{\mathbf{e}}_{y}}+B_{y,{\mathbf{y}_{f}}-2\hat{\mathbf{e}}_{y}})+\frac{9}{16}(B_{y,{\mathbf{y}_{f}}}+B_{y,{\mathbf{y}_{f}}-\hat{\mathbf{e}}_{y}})\,,\\\
\vskip 6.0pt plus 2.0pt minus 2.0pt\cr\displaystyle
B_{z,c}=-\frac{1}{16}(B_{z,{\mathbf{z}_{f}}+\hat{\mathbf{e}}_{z}}+B_{z,{\mathbf{z}_{f}}-2\hat{\mathbf{e}}_{z}})+\frac{9}{16}(B_{z,{\mathbf{z}_{f}}}+B_{z,{\mathbf{z}_{f}}-\hat{\mathbf{e}}_{z}})\,.\end{array}$
(22)

Once the pointwise values of the conservative variables are obtained by means
of Eqns. (16) and (22), the conversion to primitive variables takes place.

###  3.2 The pointwise MP5 and WENOZ spatial reconstruction algorithms

Here we introduce for the first time in a FV framework, to the extent of our
knowledge, reconstruction schemes that operate directly on the function point
values rather than on one-dimensional average quantities as traditionally done
by previous investigations (see, e.g., [40, 39, 51, 28, 31, 37]). Compared to
other FV $4^{\rm th}$-order methods, such as the one proposed by Felker &
Stone [36], we do not need to convert the average value
$\left<{U}\right>_{{\boldsymbol{c}}}$ into the mapped value
$\mathcal{V}(\left<{U_{{\boldsymbol{c}}}}\right>)$, nor to compute from it the
average value $\left<{V}\right>_{{\boldsymbol{c}}}$ before reconstructing,
thus resulting in a more efficient and cost-effective scheme. In addition, the
advantages of a point value reconstruction include that left and right states
after reconstruction can be directly fed into the Riemann solver, avoiding the
computation of interface-averaged Riemann problems. A similar feature can also
be found in the scheme presented by Núñez-de La Rosa & Munz [34], as they
reconstruct point values of the conservative quantities at the interfaces
starting from cell-averaged values via a multi-dimensional Gaussian quadrature
rule [51]. However, the approach from Núñez-de La Rosa & Munz is meant to
reconstruct either conservative or characteristic quantities (whose volume-
averaged values are the only available at the beginning of the integration
step), while in our implementation we can either reconstruct primitive or
characteristic variables (which in both approaches increase the stability of
the scheme). More importantly, the high-order reconstruction based on the
approach by [51] leads to a number of Riemann problems at each cell’s
interface equal to the number of Gaussian points used during the integration
(e.g. $4$ in the case of a three-dimensional domain when using $2$ Gaussian
points). Our scheme, on the other hand, always produces one Riemann problem
per cell’s interface.

The monotonicity preserving (MP5) scheme of Suresh & Huynh [39] provides a
$5^{\rm th}$-order accurate polynomial reconstruction at cell faces and, when
necessary, limits the reconstructed value to preserve monotonicity near
discontinuities while ensuring accuracy in smooth regions. While in the
original formulation the unlimited interface value is retrieved by one-
dimensional cell averages, our scheme replaces its definition by using point
values in the same fashion as FD schemes (see, e.g., Appendix A.2 of Del Zanna
et al. [20] for a brief summary). The unlimited value at $i+\frac{1}{2}$ is
now given by

$\left.\ P(x)\right|_{i+\frac{1}{2}}\equiv
V^{L}_{i+\frac{1}{2}}=\frac{1}{128}\left(3V_{i-2}-20V_{i-1}+90V_{i}+60V_{i+1}-5V_{i+2}\right)\,.\vspace{0.2
cm}$ (23)

where $\\{V_{i\pm 2},V_{i\pm 1},V_{i}\\}$ are now primitive quantities
evaluated at cell-centers. The scheme employs, as in the original case, a
five-point stencil in order to distinguish between local extrema and a genuine
O(1) discontinuity. The MP5 limiting procedure remains unchanged as
monotonicity-preserving bounds are equally well expressed in terms of point
values. We refer the reader to the original work of Suresh & Huynh [39] or to
Appendix 2 of Del Zanna et al. [20].

Likewise, we also introduce a pointwise version of the $5^{\rm th}$-order
weighted essentially non-oscillatory (WENO, [40]) scheme in the version
proposed by Borges et al. (WENOZ, [41]), which attains better spatial
resolution with reduced numerical dissipation at a modest computational cost.
The same interface value of Eq. (23) at $x=x_{i+\frac{1}{2}}$ is computed as
the convex combination of $3^{\rm rd}$-order accurate interface values built
on the three possible three-point sub-stencils $\\{i-2,i-1,i\\}$,
$\\{i-1,i,i-1\\}$, and $\\{i,i+1,i+2\\}$, i.e.

$\displaystyle
V^{L}_{i+\frac{1}{2}}=\displaystyle\omega_{0}\frac{3V_{i-2}-10V_{i-1}+15V_{i}}{8}+\omega_{1}\frac{-V_{i-1}+6V_{i}+3V_{i+1}}{8}+\omega_{2}\frac{3V_{i}+6V_{i+1}-V_{i+2}}{8}\,,$
(24)

where, again, $\\{V_{i\pm 2},V_{i\pm 1},V_{i}\\}$ are point values. The
weights $\omega_{l}$, for $l=\\{0,1,2\\}$, are defined as

$\omega_{l}=\frac{\alpha_{l}}{\sum_{m}\alpha_{m}},\hskip
14.22636pt\alpha_{l}=d_{l}\left(1+\frac{|\beta_{0}-\beta_{2}|}{\beta_{l}+\epsilon}\right)\,.$
(25)

Here $d_{l}$ denotes the optimal weights that reckon the $5^{\rm th}$-order
accurate approximation of Eq. (23), while $\epsilon=10^{-40}$ is a small
number which avoids division by zero. In our new framework, the optimal
weights $\\{d_{0}=1/16,d_{1}=5/8,d_{2}=5/16\\}$ can be determined from the
uniqueness of the interpolating polynomial with $d_{0}+d_{1}+d_{2}=1$. On the
other hand, the smoothness indicators $\beta_{l}$ estimate the regularity of
the polynomial approximation and, for a reconstruction based on point values,
have been found to be the same as in Borges [41]:

$\displaystyle\beta_{0}$
$\displaystyle=\frac{13}{12}(V_{i-2}-2V_{i-1}+V_{i})^{2}+\frac{1}{4}(V_{i-2}-4V_{i-1}+3V_{i})^{2}\,,$
(26) $\displaystyle\beta_{1}$
$\displaystyle=\frac{13}{12}(V_{i-1}-2V_{i}+V_{i+1})^{2}+\frac{1}{4}(V_{i-1}-V_{i+1})^{2}\,,$
$\displaystyle\beta_{2}$
$\displaystyle=\frac{13}{12}(V_{i}-2V_{i+1}+V_{i+2})^{2}+\frac{1}{4}(3V_{i}-4V_{i+1}+V_{i+2})^{2}\,.$

Note that the right state, $V_{i-\frac{1}{2}}^{R}$, is obtained similarly by
reversing the index order (e.g. $i-1\to i+1$, etc).

We will also occasionally employ the $3^{\rm rd}$-order WENO scheme (here
denoted as WENO3) [40] in the improved version of Yamaleev & Carpenter [52].
The point value approach leads to the following expression for the left
interface value

$V_{i+\frac{1}{2}}^{L}=\omega_{0}\frac{V_{i+1}+V_{i}}{2}+\omega_{1}\frac{-V_{i+1}+3V_{i}}{2}\,,$
(27)

where the weights are given by:

$\omega_{l}=\frac{\alpha_{l}}{\alpha_{0}+\alpha_{1}}\,,\quad\alpha_{l}=d_{l}\left(1+\frac{|\Delta_{i+\frac{1}{2}}-\Delta_{i-\frac{1}{2}}|^{2}}{\beta_{l}+\epsilon}\right)\,,\quad\beta_{0}=\Delta_{i-\frac{1}{2}}^{2}\,,\quad\beta_{1}=\Delta_{i+\frac{1}{2}}^{2}\,,$
(28)

with $l=\\{0,1\\}$. In Eq. (28), $\Delta_{i+\frac{1}{2}}=(V_{i+1}-V_{i})$
whereas $\\{d_{0}=3/4$, $d_{1}=1/4\\}$ (note that in the standard finite
volume formulation, starting from 1D cell averages, one has $d_{0}=2/3$,
$d_{1}=1/3$). Finally, in order to avoid loss of accuracy at critical points
(see [52] and [15]) we adopt here $\epsilon=\Delta x^{2}$.

###  3.3 Flux computation and temporal integration

Spatial reconstruction provides the $4^{\rm th}$-order accurate pointwise
values for left/right states at a given face center (e.g.
$V^{L}_{{\mathbf{x}_{f}}}$,$V^{R}_{{\mathbf{x}_{f}}}$). These serve as input
states to the Riemann solver for the computation of the high-order interface
fluxes. Our implementation comes with a variety of Riemann solvers, such as
the linearized Roe solver [53], HLL-type Riemann solvers [54], including the
HLLC [55, 56, 11] and HLLD [57, 58], the Lax-Friedrichs [59], and the GFORCE
[60] solver. Denoting with $\mathcal{F}_{x,y,z}(\cdot,\,\cdot)$ a generic 1-D
Riemann solver along the $\\{x,y,z\\}$ direction, the interface upwind
punctual fluxes are obtained as:

$F_{x,{\mathbf{x}_{f}}}=\mathcal{F}_{x}(V^{L}_{\mathbf{x}_{f}},V^{R}_{\mathbf{x}_{f}})\,,\quad
F_{y,{\mathbf{x}_{f}}}=\mathcal{F}_{y}(V^{L}_{\mathbf{y}_{f}},V^{R}_{\mathbf{y}_{f}})\,,\quad
F_{z,{\mathbf{x}_{f}}}=\mathcal{F}_{z}(V^{L}_{\mathbf{z}_{f}},V^{R}_{\mathbf{z}_{f}})\,.$
(29)

$4^{\rm th}$-order accurate surface averages are then obtained by the
quadrature rule in Eq. (17) using the transverse Laplacian operators:

$\displaystyle\hat{F}_{x,{\mathbf{x}_{f}}}=F_{x,{\mathbf{x}_{f}}}+\frac{\Delta^{x}_{\bot}F_{x,{\mathbf{x}_{f}}}}{24},\quad\hat{F}_{y,{\mathbf{y}_{f}}}=F_{y,{\mathbf{x}_{f}}}+\frac{\Delta^{y}_{\bot}F_{y,{\mathbf{x}_{f}}}}{24},\quad\hat{F}_{z,{\mathbf{z}_{f}}}=F_{z,{\mathbf{x}_{f}}}+\frac{\Delta^{z}_{\bot}F_{z,{\mathbf{x}_{f}}}}{24}\,.$
(30)

Note that this requires an additional integration in the first layer of ghost
zones in order to retrieve the flux information. Moreover, thanks to the
deployment of point value reconstructions, we do not need to solve a Riemann
problem with interface-averaged primitive states as it is instead done by
Felker & Stone (Eq. 17 in [36]).

These fluxes are used to build the right-hand side of Eq. (10) which will be
hereafter denoted with $\mathcal{L}(U)$. Our numerical scheme employs the five
stage $4^{\rm th}$-order explicit Strong Stability Preserving Runge-Kutta
method (eSSPRK(5,4), [43, 44]) to discretize Eq. (10)

$\begin{array}[]{lcl}U^{(1)}&=&\displaystyle U^{n}+0.391752226571890\Delta
t\mathcal{L}(U^{n})\\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr
U^{(2)}&=&\displaystyle
0.444370493651235U^{n}+0.555629506348765U^{(1)}+0.368410593050371\Delta
t\mathcal{L}(U^{(1)})\\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr
U^{(3)}&=&\displaystyle
0.620101851488403U^{n}+0.379898148511597U^{(2)}+0.251891774271694\Delta
t\mathcal{L}(U^{(2)})\\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr
U^{(4)}&=&\displaystyle
0.178079954393132U^{n}+0.821920045606868U^{(3)}+0.544974750228521\Delta
t\mathcal{L}(U^{(3)})\\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr
U^{n+1}&=&\displaystyle
0.517231671970585U^{(2)}+0.096059710526147U^{(3)}+0.063692468666290\Delta
t\mathcal{L}(U^{(3)})\\\ \vskip 6.0pt plus 2.0pt minus
2.0pt\cr&&\displaystyle+0.386708617503268U^{(4)}+0.226007483236906\Delta
t\mathcal{L}(U^{(4)})\,.\end{array}$ (31)

The eSSPRK(5,4) method is strongly stable (namely, no bounded temporal growth
is allowed: $\|U^{n+1}\|\leq\|U^{n}\|$) under the assumption that the forward
Euler method employed is strongly stable under the Courant-Friedrichs-Levy
(CFL) restriction ($\|U^{n}+\Delta t\mathcal{L}(U^{n})\|\leq\|U^{n}\|$, [61,
62]). In our implementation, the CFL condition determines the time step
according to

$\Delta t={\cal
C}\displaystyle\min_{ijk}\left[\displaystyle\min_{d}\left(\displaystyle\frac{\Delta
x^{d}}{|\lambda^{d}_{\rm max}|}\right)\right]$ (32)

where ${\cal C}$ is the CFL number, $\lambda^{d}_{\rm max}$ is the outermost
fast magnetotosonic wave and $\Delta x^{d}$ is the cell length through the
different directions $d$.

For a detailed description of the stability properties of the method we refer
the reader to [62, 44, 63, 64, 65].

###  3.4 Order reduction in the presence of discontinuities

The recovery of point values from volume averages follows the assumption that
the function $Q$ is locally continuous up to the $3^{\rm rd}$ derivative in
its expansion. However, this assumption manifestly breaks down in the
immediate neighbour of a discontinuity, where using Eq. (16) may easily
produce unphysical values such as negative energies or densities. In order to
address this shortcoming, we modify the recovery of pointwise conservative
variables at the beginning of any integration stage by modifying Eq. (16) as

$U_{\boldsymbol{c}}=\displaystyle\left<{U}\right>_{\boldsymbol{c}}-\theta_{\boldsymbol{c}}\frac{\Delta\left<{U}\right>_{\boldsymbol{c}}}{24}+O(h^{4})\,,$
(33)

where $\theta_{\boldsymbol{c}}$ has been introduced to detect the presence of
a discontinuity

$\theta_{\boldsymbol{c}}=\left\\{\begin{array}[]{ll}1&\quad{\rm
if}\quad\eta_{\boldsymbol{c}}<\eta_{d}\,,\\\ \vskip 6.0pt plus 2.0pt minus
2.0pt\cr 0&\quad{\rm otherwise}\,.\end{array}\right.$ (34)

The parameter $\eta_{\boldsymbol{c}}$ is defined as
$\eta_{\boldsymbol{c}}=\sqrt{\eta^{2}_{{\boldsymbol{c}},x}+\eta^{2}_{{\boldsymbol{c}},y}+\eta^{2}_{{\boldsymbol{c}},z}}$,
while $\eta_{d}$ is a threshold parameter. When $\theta_{\boldsymbol{c}}=0$ in
Eq. (33), the local order of the scheme is reduced in proximity of
discontinuous fronts, including avoiding the integration of fluxes (Eq. 30).
In these regions, we further increase the robustness of the scheme by also
lowering the order of the reconstruction scheme by switching to the $3^{\rm
rd}$-order WENO scheme (see end of §3.2) or to piecewise linear
reconstruction. This approach (known as the “fallback”) has been employed, in
a similar fashion, also by [34] in the context of high-order schemes. We now
propose two different discontinuity detectors in order to evaluate
$\eta_{\boldsymbol{c}}$.

The first one is based on Jameson’s pressure-based shock sensor [66]

$\eta_{{\boldsymbol{c}},x}=\frac{\left<{Q}\right>_{{\boldsymbol{c}}+\hat{\mathbf{e}}_{x}}-2\left<{Q}\right>_{{\boldsymbol{c}}}+\left<{Q}\right>_{{\boldsymbol{c}}-\hat{\mathbf{e}}_{x}}}{|\left<{Q}\right>_{{\boldsymbol{c}}+\hat{\mathbf{e}}_{x}}|+2|\left<{Q}\right>_{{\boldsymbol{c}}}|+|\left<{Q}\right>_{{\boldsymbol{c}}-\hat{\mathbf{e}}_{x}}|+\epsilon}\,,$
(35)

(similar expressions hold for $\eta_{{\boldsymbol{c}},y}$ and
$\eta_{{\boldsymbol{c}},z}$) where $\epsilon$ is a small number to prevent
negligible variations from triggering the switch. Eq. (35) is first evaluated
with $\left<{Q}\right>$ being the volume-averaged density, pressure, and
magnetic pressure and the maximum over these three values is finally
considered. Jameson’s switch has the advantage of keeping the stencil compact
(no additional ghost zones are needed), albeit three zones may not be
sufficient to distinguish a discontinuity from a smooth extremum [39].

The second discontinuity sensor that we propose employs a derivatives ratio
that attempts to identify high-frequency oscillations which typically arise
near a discontinuous front. Since high wave-number Fourier modes will not be
dumped around these critical points, we expect the contribution of higher
order derivatives to be non-negligible. To this end, we first evaluate the
normalized contribution of odd and even modes separately by evaluating the
parameters $\eta^{o}_{\boldsymbol{c}}$ and $\eta^{e}_{\boldsymbol{c}}$

$\eta^{o}_{\boldsymbol{c}}=\frac{|\delta^{(3)}\left<{Q}\right>_{\boldsymbol{c}}|}{|\left<{Q}\right>_{{\boldsymbol{c}},{\rm
ref}}|+|\delta^{(1)}\left<{Q}\right>_{\boldsymbol{c}}|+|\delta^{(3)}\left<{Q}\right>_{\boldsymbol{c}}|+\epsilon}\,,\qquad\eta^{e}_{\boldsymbol{c}}=\frac{|\delta^{(4)}\left<{Q}\right>_{\boldsymbol{c}}|}{|\left<{Q}\right>_{{\boldsymbol{c}},{\rm
ref}}|+|\delta^{(2)}\left<{Q}\right>_{\boldsymbol{c}}|+|\delta^{(4)}\left<{Q}\right>_{\boldsymbol{c}}|+\epsilon}\,,$
(36)

where $\delta^{(m)}Q_{\boldsymbol{c}}$ (with $m=1,..,4$) are undivided $2^{\rm
nd}$-order accurate approximations to the $m$-th derivative. In the
$x$-direction, for instance:

$\begin{array}[]{lcl}\delta^{(1)}\left<{Q}\right>_{\boldsymbol{c}}&=&\displaystyle\frac{1}{2}\left(\left<{Q}\right>_{{\boldsymbol{c}}+\hat{\mathbf{e}}_{x}}-\left<{Q}\right>_{{\boldsymbol{c}}-\hat{\mathbf{e}}_{x}}\right)\,,\\\
\vskip 6.0pt plus 2.0pt minus
2.0pt\cr\delta^{(2)}\left<{Q}\right>_{\boldsymbol{c}}&=&\displaystyle{\left<{Q}\right>_{{\boldsymbol{c}}+\hat{\mathbf{e}}_{x}}-2\left<{Q}\right>_{\boldsymbol{c}}+\left<{Q}\right>_{{\boldsymbol{c}}-\hat{\mathbf{e}}_{x}}}\,,\\\
\vskip 6.0pt plus 2.0pt minus
2.0pt\cr\delta^{(3)}\left<{Q}\right>_{\boldsymbol{c}}&=&\displaystyle\frac{1}{2}\left(\left<{Q}\right>_{i+2\hat{\mathbf{e}}_{x}}-2(\left<{Q}\right>_{{\boldsymbol{c}}+\hat{\mathbf{e}}_{x}}-2\left<{Q}\right>_{{\boldsymbol{c}}-\hat{\mathbf{e}}_{x}})-\left<{Q}\right>_{{\boldsymbol{c}}-2\hat{\mathbf{e}}_{x}}\right)\,,\\\
\vskip 6.0pt plus 2.0pt minus
2.0pt\cr\delta^{(4)}\left<{Q}\right>_{\boldsymbol{c}}&=&\displaystyle\left<{Q}\right>_{{\boldsymbol{c}}+2\hat{\mathbf{e}}_{x}}-4\left<{Q}\right>_{{\boldsymbol{c}}+\hat{\mathbf{e}}_{x}}+6\left<{Q}\right>_{\boldsymbol{c}}-4\left<{Q}\right>_{{\boldsymbol{c}}-\hat{\mathbf{e}}_{x}}+\left<{Q}\right>_{{\boldsymbol{c}}-2\hat{\mathbf{e}}_{x}}\,.\end{array}$
(37)

In our experience, using primitive variables has shown to yield a more robust
identification method. For this reason, we take
$\left<{Q}\right>_{{\boldsymbol{c}}}={\cal
V}(\left<{U}\right>_{\boldsymbol{c}})$ at the beginning of the Runge-Kutta
stage. The reference value $\left<{Q}\right>_{{\boldsymbol{c}},\rm ref}$ in
Eq. (36) is equal to $\left<{Q}\right>_{{\boldsymbol{c}}}$ for density and
pressure while we employ $\sqrt{p_{\boldsymbol{c}}/\rho_{\boldsymbol{c}}}$
(isothermal sound speed) for velocity components and $|B|_{\boldsymbol{c}}$
for magnetic field components. Once $\eta^{o}_{\boldsymbol{c}}$ and
$\eta^{e}_{\boldsymbol{c}}$ have been computed, we simply take the maximum of
the two:
$\eta_{\boldsymbol{c}}=\max(\eta_{\boldsymbol{c}}^{o},\eta_{\boldsymbol{c}}^{e})$.

The derivative ratio detector needs a larger stencil thereby asking for one
additional ghost zone. Along the same arguments, a more sophisticated approach
that requires an even larger stencil, is shown by Bambozzi & Pires [67].

##  4 The high-order constrained transport

Figure 1: Left: positioning of MHD variables in the UCT formalism. Staggered
magnetic field components (blue) are face-centered, the electromotive force is
edge-centered (red) while remaining HD quantities are located at the zone
center (green). Right: top view of the intersection between four neighbor
zones: N, S, E, and W indicate the four cardinal directions with respect to
the zone edge (here represented by the intersection between four neighbor
zones), $\altmathcal{R}_{x}(F_{{\mathbf{y}_{f}}})$ and
$\altmathcal{R}_{y}(F_{{\mathbf{x}_{f}}})$ are 1-D reconstruction operators
applied to each zone face, see Eq. (38).

As we outlined in §2, the HD sub-system cannot be straightforwardly extended
to the MHD system without introducing numerical strategies to control the
rising of spurious magnetic monopoles. Failure to fulfill this requisite may
lead to unphysical effects such as plasma acceleration parallel to the field,
incorrect jump conditions, wrong propagation speed of discontinuities and odd-
even decoupling [48, 6, 5]. In what follows we employ the upwind constrained
transport method [45, 3, 5, 46, 47, 18] in order to ensure that the
divergence-free condition is maintained to machine accuracy.

In the FV-CT formalism the field’s surface integrals
$\hat{B}_{x,{\mathbf{x}_{f}}}$, $\hat{B}_{y,{\mathbf{y}_{f}}}$,
$\hat{B}_{z,{\mathbf{z}_{f}}}$ are the primary representation of the magnetic
field evolved in time. Eq. (4) regulates the time variation of the magnetic
field fluxes, implying that the formal correspondence

$\begin{array}[]{lll}F^{[B_{x}]}_{x}=0\,,&F^{[B_{x}]}_{y}=E_{z}\,,&F^{[B_{x}]}_{z}=-E_{y}\,,\\\
\vskip 6.0pt plus 2.0pt minus 2.0pt\cr
F^{[B_{y}]}_{x}=-E_{z}\,,&F^{[B_{y}]}_{y}=0\,,&F^{[B_{y}]}_{z}=E_{x}\,,\\\
\vskip 6.0pt plus 2.0pt minus 2.0pt\cr
F^{[B_{z}]}_{x}=E_{y}\,,&F^{[B_{z}]}_{y}=-E_{x}\,,&F^{[B_{z}]}_{z}=0\,,\end{array}$
(38)

holds when solving a Riemann problem. While this certainly holds for the
fluxes computed at the face midpoint, Eq. (14) requires the EMF to become
available at a zone edge to serve as a discrete version of Stokes theorem. In
order to reconstruct the induction fluxes from the face center (where they are
available as point values) to one of the four adjacent edges, we employ the
UCT method of Mignone & Del Zanna [18]. This method provides a general
formalism for systematic construction of EMF averaging procedures where
dissipation terms depend on the left and right transverse magnetic field
components.

Referring to the right panel of Fig. 1 and being consistent with the notation
presented in Mignone & Del Zanna [18], we consider a top view of the
intersection of four zones at a cell-edge. We denote, respectively, the left
and right reconstructed states along the y-coordinate as south (S) and north
(N), while on the x-direction the reconstructed states are labeled with west
(W) or east (E) with respect to the intersection point. The edge-centered EMF
with the desired upwind properties can be build up from

$E_{z,{\mathbf{z}_{e}}}=-\left[\left(a_{x}\overline{v}_{x}B_{y}\right)^{W}+\left(a_{x}\overline{v}_{x}B_{y}\right)^{E}\right]+\left[\left(a_{y}\overline{v}_{y}B_{x}\right)^{N}+\left(a_{y}\overline{v}_{y}B_{x}\right)^{S}\right]+\left[\left(d_{x}B_{y}\right)^{E}-\left(d_{x}B_{y}\right)^{W}\right]-\left[\left(d_{y}B_{x}\right)^{N}-\left(d_{y}B_{x}\right)^{S}\right]\,,$
(39)

where $\overline{v}_{x,y}$ are the transverse velocities reconstructed from
the interfaces, $a_{x,y}$ and $d_{x,y}$ are the flux and diffusion
coefficients, and $B_{x,y}$ are the pointwise values of the magnetic field
reconstructed from the corresponding interface. The other EMF components may
be obtained by cyclic permutations. For the sake of briefness, we address the
reader to find all the details of the different EMF averaging techniques that
we employ in the present work in Mignone & Del Zanna [18]. Note that the
employment of point value reconstruction methods (rather than 1D volume
average-based ones) requires the electric field to be evaluated over 3 layers
of ghost zones in the transverse direction rather than 4.

Finally, the line-averaged integrated electric field is computed as

$\bar{E}_{x,{\mathbf{x}_{e}}}=\displaystyle
E_{x,{\mathbf{x}_{e}}}+\frac{\Delta^{x}E_{x,{\mathbf{x}_{e}}}}{24},\quad\bar{E}_{y,{\mathbf{y}_{e}}}=\displaystyle
E_{y,{\mathbf{y}_{e}}}+\frac{\Delta^{y}E_{y,{\mathbf{y}_{e}}}}{24},\quad\bar{E}_{z,{\mathbf{z}_{e}}}=\displaystyle
E_{z,{\mathbf{z}_{e}}}+\frac{\Delta^{z}E_{z,{\mathbf{z}_{e}}}}{24}\,,$ (40)

and then employed to update the staggered magnetic field components using the
discrete version of Stoke’s theorem Eq. (14).

In the following we recap the steps of our $4^{\rm th}$-order numerical
scheme.

  1. 1.

Start with the volume average of conserved quantities
$\left<{U}\right>_{{\boldsymbol{c}}}$ and the face average of the staggered
magnetic field components, $\hat{B}_{x,{\mathbf{x}_{f}}}^{x}$,
$\hat{B}_{y,{\mathbf{y}_{f}}}^{y}$, $\hat{B}_{z,{\mathbf{z}_{f}}}^{z}$.

  2. 2.

Assign boundary conditions to these variables.

  3. 3.

Obtain volume-averaged primitive variables by converting the volume averages
$\left<{V}\right>_{c}={\cal V}(\left<{U}\right>_{\boldsymbol{c}})$. Notice
that this array is $2^{\rm nd}$-order accurate and it is only used in the
limiting process, thus it does not replace the pointwise conversion later on.

  4. 4.

Evaluate the $\theta_{\boldsymbol{c}}$ parameter for order reduction in
presence of discontinuities by means of either the Jameson’s shock sensor (Eq.
35) or by weighing up the contribution of higher-order derivatives (Eq. 36)
and flag troubled cells along with first neighbours.

  5. 5.

Get pointwise value $U_{\boldsymbol{c}}$ (Eq. 33) at the cell center and
$B_{{\mathbf{x}_{f}}},B_{{\mathbf{y}_{f}}},B_{{\mathbf{z}_{f}}}$ (Eq. 21) at
the corresponding face centers. Also, during this step, interpolate the
magnetic field from the face center to the cell center (Eq. 22).

  6. 6.

Convert pointwise conservative variables $U_{\boldsymbol{c}}$ in primitive
variables $V_{\boldsymbol{c}}={\cal V}(U_{\boldsymbol{c}})$.

  7. 7.

Reconstruct L/R interface values (Eq. 23 or 24) in primitive (or
characteristic) variables.

  8. 8.

Obtain upwind fluxes of hydro quantities (Eq. 29) and average them over the
corresponding faces (Eq. 30). Concurrently, store the transverse velocities
(e.g., $\overline{v}_{y}$ and $\overline{v}_{z}$ at an $x$-interface) as well
as the coefficients $a^{L/R}$ and $d^{L/R}$ available with the chosen Riemann
solver at a zone interface. These quantities are needed in the next step of
the UCT scheme. For completeness we report the explicit form of these
coefficients in B.

  9. 9.

Reconstruct the transverse velocities obtained with the 1D Riemann solver and
the normal components of the magnetic field from the face centers to the
edges. This will yield the $N,S$ and $W,E$ values in Eq. (39). At an
$x$-interface, for instance, one has to reconstruct
$\overline{v}_{y,{\mathbf{x}_{f}}}$ and $B_{x,{\mathbf{x}_{f}}}$ along the
$y$-direction, and then, again, $\overline{v}_{z,{\mathbf{x}_{f}}}$ and
$B_{x,{\mathbf{x}_{f}}}$ along the $z$-direction.

  10. 10.

Obtain the electric field at corners (Eq. 39) by means of different averaging
procedures (see B).

  11. 11.

Line-average the electric field according to Eq. (40).

  12. 12.

Add all contributions to the right-hand side and obtain the solution
$\left<{U}\right>_{{\boldsymbol{c}}}$ at the next time level or RK stage (Eq.
31).

##  5 Numerical benchmarks

We now test the capabilities of the novel method in a series of 2D and 3D
numerical benchmarks. Each test is validated twice, by solving for the MHD and
the RMHD case. As a reference, we will compare our results against the
traditional $2^{\rm nd}$\- or $3^{\rm rd}$-order schemes. As a default, Eq.
(10) will be discretized by the $2^{\rm nd}$-, $3^{\rm rd}$-, or $4^{\rm
th}$-order schemes by means of a RK2, RK3, or RK4 time integrator. The WENOZ
and MP5 spatial reconstructions will always use point values and, otherwise
differently stated, each MHD and RMHD benchmark will employ, respectively, a
HLLD solver coupled with the UCT-HLLD EMF average, and the GFORCE solver plus
the UCT-GFORCE average that we suitably extended to relativistic frameworks.

Whenever the initial condition is a smooth analytic function our scheme
integrates it by means of a Gaussian quadrature rule employing 4 Gaussian
points.

We evaluate $L_{1}$ norm errors for a generic quantity $Q$ against a reference
solution $Q_{\rm ref}$ as $\epsilon_{1}(Q)=\sum_{i}|Q(x_{i})-Q_{\rm
ref}(x_{i})|/N$, where $x_{i}$ is a generic point of the domain and $N$ is the
total number of grid zones. The convergence rate $O_{L_{1}}$ is thus
calculated via a linear fit of $\log(\epsilon_{1})$ against $\log(N_{x})$,
with $N_{x}$ being the linear resolution.

All computations have been carried out by means of the PLUTO astrophysical
code [68, 69], where the $4^{\rm th}$-order method has been implemented.

###  5.1 3D Rotated Shock Tube

Figure 2: Primitive variable profiles for the 3D rotated MHD shock tube
problem at $t_{f}=t^{\prime}_{f}\cos\alpha\cos\gamma$, where
$t^{\prime}_{f}=0.2$ obtained with the WENOZ $4^{\rm th}$-order scheme (red
squares) on a uniform grid of $768\times 8\times 8$ zones. All the variables
are plotted along the non-rotated direction $x^{\prime}$. On the left panels
density (top) and pressure (bottom) are shown, while in the remaining panels,
the three components of velocity (top) and magnetic field (bottom) are
plotted. The black solid lines represent a reference 1D solution computed by
emplyoing a $2^{\rm nd}$-order scheme with a resolution of $N_{x}=16384$ grid
points. Figure 3: Same of Fig. 2 but for the relativistic case. Here the
numerical scheme adopted is using the MP5 algorithm and the “unrotated” final
time is $t_{f}^{\prime}=0.55$.

We begin this section by considering a 3D rotated version of a standard shock-
tube problem. The initial condition is best expressed in the original (non-
rotated) frame where $x^{\prime}$ is the coordinate perpendicular to the plane
of discontinuity. In the MHD case we have

$\left\\{\begin{array}[]{lclr}\altmathcal{V}^{\prime}_{L}&=&\displaystyle\Big{(}1.08,\,1.2,\,0.01,\,0.5,\,\frac{2}{\sqrt{4\pi}},\,\frac{3.6}{\sqrt{4\pi}},\,\frac{2}{\sqrt{4\pi}},\,0.95\Big{)}&\quad{\rm
for}\quad x^{\prime}<0\,,\\\ \vskip 6.0pt plus 2.0pt minus
2.0pt\cr\altmathcal{V}^{\prime}_{R}&=&\displaystyle\Big{(}1,\,0,\,0,\,0,\,\frac{2}{\sqrt{4\pi}},\,\frac{4}{\sqrt{4\pi}},\,\frac{2}{\sqrt{4\pi}},\,1\Big{)}&\quad{\rm
for}\quad x^{\prime}>0\,,\end{array}\right.$ (41)

where
$\altmathcal{V}^{\prime}=(\rho,v^{\prime}_{x^{\prime}},v^{\prime}_{y^{\prime}},v^{\prime}_{z^{\prime}},B^{\prime}_{x^{\prime}},B^{\prime}_{y^{\prime}},B^{\prime}_{z^{\prime}},p)$
is an array of primitive quantities in the non-rotated frame. In the
relativistic case we have, likewise,

$\left\\{\begin{array}[]{lclr}\altmathcal{V}^{\prime}_{L}&=&\displaystyle\left(1.08,\,0.4,\,0.3,\,0.2,\,2,\,0.3,\,0.3,\,0.95\right)&\quad{\rm
for}\quad x^{\prime}<0\,,\\\ \vskip 6.0pt plus 2.0pt minus
2.0pt\cr\altmathcal{V}^{\prime}_{R}&=&\displaystyle\left(1,\,-0.45,\,-0.2,\,0.2,\,2,\,-0.7,\,0.5,\,1\right)&{\rm
for}\quad x^{\prime}>0\,.\end{array}\right.$ (42)

For each case, we rotate the initial condition by an angle $\alpha$ around the
$z$-axis and then by an angle $\gamma$ around the $y$-axis. The coordinate
transformation used for the 3D rotation is explained in detail in [15] (see
§4.2.2 in that paper) and leads to a planar symmetry condition whereby, for
any quantity $q(\mathbf{x})$, one must have
$q(\mathbf{x}+\mathbf{s})=q(\mathbf{x})$, where $\mathbf{s}=(n_{x}\Delta
x,\,n_{y}\Delta y,\,n_{z}\Delta z)$ corresponds to an integer shift of cells.
We choose $\tan\alpha=1/2$ and $\tan\gamma=\cos\alpha\tan\beta=-1/(2\sqrt{5})$
with $\tan\beta=-1/4$ corresponding to a translational invariance by an
integer shift of cells, i.e., $(n_{x},\,n_{y},\,n_{z})=(1,-2,0)$ and
$(n_{x},\,n_{y},\,n_{z})=(1,0,4)$.

We solve the MHD / RMHD equations using our $4^{\rm th}$-order method,
adopting the WENOZ reconstruction in characteristic variables for the
classical version and the MP5 reconstruction in primitive variables for the
relativistic variant. The computational domain is the box
$[-L/2,\,L/2]\times[0,\,L/96]\times[0,\,L/96]$ covered by a uniform grid of
$768\times 8\times 8$ zones. We evolve the solution until
$t_{f}=t^{\prime}_{f}\cos\alpha\cos\gamma=4/\sqrt{21}$ where
$t^{\prime}_{f}=0.2$ (for classical MHD), $t^{\prime}_{f}=0.55$ is the final
time for the non-rotated (1D) versions of the problem, see, e.g., [70, 71, 15,
15] for the classical MHD version and [72, 58, 73] (and references therein)
for the relativistic variant. One dimensional plots of relevant fluid
quantities are shown in the unrotated frame in Fig. (2) and (3) against a 1D
high-resolution reference profile computed with a $2^{\rm nd}$-order scheme
with a resolution of $N_{x}=16384$ grid points. The solution features a
contact discontinuity in the middle ($x\approx 0$) separating two outermost
fast magnetosonic shocks enclosing a pair of rotational waves followed by two
slow waves. Profiles are correctly reproduced and oscillations in the normal
component of magnetic field are of the same order as those of [15] and smaller
than those in [71] who used a $2^{\rm nd}$-order scheme. Our results
demonstrate that our $4^{\rm th}$-order method can correctly capture
propagating discontinuous fronts although the order of the scheme degenerates
to $1^{\rm st}$-order.

###  5.2 2D iso-density and iso-enthalpy vortex

Figure 4: Zoomed plot of the thermal pressure profile in the 2D MHD iso-
density vortex (left panel) and density in the 2D RMHD iso-enthalpy vortex
(right panel). In each panel the exact solution’s profile (continuous black
line) is compared with the profiles evolved by a $4^{\rm th}$-order scheme
adopting the WENOZ (red squares), the MP5 (blue circles), and the $2^{\rm
nd}$-order LINEAR scheme (green triangles). Filled symbols refer to data
produced at a resolution of $N_{x}=128$, while the empty ones are data at
$N_{x}=512$.

We consider the MHD iso-density vortex advection problem (as introduced by
Balsara [7] and later presented by Mignone [15] and Dumbser [74]) and also
propose a new exact equilibrium solution of the ideal RMHD equations with
constant enthalpy. For the MHD case, the initial condition consists of a
magnetized vortex at constant density in force balance that propagates along
the main diagonal of a computational domain of $[-5,5]\times[-5,5]$ under
periodic boundary conditions. The equilibrium condition is best expressed in
cylindrical coordinates $(R,\phi,z)$, from radial momentum balance

$\displaystyle\frac{d}{dR}\left(p+\frac{B_{z}^{2}}{2}\right)=\frac{\rho
v^{2}_{\phi}}{R}-\frac{1}{2R^{2}}\frac{d}{dR}\left(R^{2}B^{2}_{\phi}\right)\,.$
(43)

We assume constant density $\rho=1$ and Eq. (43) is satisfied by adopting the
following profiles:

$\mathbf{v}=\mathbf{v}_{0}+(-y,x)ke^{q(1-R^{2})}\,,\quad\mathbf{B}=(-y,x)\mu
e^{q(1-R^{2})}+B_{z}\hat{\mathbf{e}}_{z}\,,$ (44)

where $\mathbf{v}_{0}$ is the translational velocity, $R=\sqrt{x^{2}+y^{2}}$
is the cylindrical radius, $k=\mu=1/2\pi$ and the parameter $q$ is set to $1$
as in [74, 15]. The pressure term on the left hand side of Eq. (43) is then
readily obtained from

$p_{\rm t}\equiv
p+\frac{B_{z}^{2}}{2}=1+\frac{1}{4q}\left[\mu^{2}\left(1-2qR^{2}\right)-fk^{2}\right]e^{2q(1-R^{2})}\,,$
(45)

where $f=\rho$ for the classical case. We equally distribute $p_{t}$ in Eq.
(45) among the thermal and magnetic contributions, i.e. $p=\alpha p_{\rm t}$
and $B_{z}^{2}/2=(1-\alpha)p_{\rm t}$, where $\alpha\in[0,1]$. In this setup,
we use $\alpha=1/2$ and $\mathbf{v}_{0}=(1,1)$.

In the RMHD case (for which we also assume a purely radial dependence), the
time-independent radial component of the momentum equation takes instead the
form

$\displaystyle\frac{dp}{dR}=\frac{w\gamma
v_{\phi}^{2}}{R}+\frac{E_{R}}{R}\frac{d\left(RE_{R}\right)}{dR}-B_{z}\frac{dB_{z}}{dR}-\frac{B_{\phi}}{R}\frac{d\left(RB_{\phi}\right)}{dR}\,,$
(46)

where $w=5$ is the constant enthalpy. We assume $B_{z}=E_{R}=0$, so that Eq.
(46) has solution given by Eq. (44) and (45) with $f=w$. Here we choose a
static vortex ($v_{0}=0$). Density is recovered from the ideal equation of
state (adiabatic index $\Gamma=5/3$) as $\rho=\displaystyle w-5p/2$, which is
ensured to be strictly positive by our choice of the initial enthalpy.

$N_{x}$ | MHD MP5 |  | MHD WENOZ |  | RMHD MP5 |  | RMHD WENOZ
---|---|---|---|---|---|---|---
| $\epsilon_{1}(B_{x})$ | $O_{L_{1}}$ |  | $\epsilon_{1}(B_{x})$ | $O_{L_{1}}$ |  | $\epsilon_{1}(B_{x})$ | $O_{L_{1}}$ |  | $\epsilon_{1}(B_{x})$ | $O_{L_{1}}$
$32$ | 9.99E-04 | - |  | 9.01E-04 | - |  | 3.93E-04 | - |  | 3.11E-04 | -
$64$ | 6.09E-05 | 4.04 |  | 6.00E-05 | 3.91 |  | 1.70E-05 | 4.53 |  | 1.83E-05 | 4.09
$128$ | 2.35E-06 | 4.70 |  | 2.37E-06 | 4.66 |  | 6.09E-07 | 4.80 |  | 6.41E-07 | 4.84
$256$ | 1.06E-07 | 4.47 |  | 1.06E-07 | 4.48 |  | 2.34E-08 | 4.71 |  | 2.36E-08 | 4.77
$512$ | 5.76E-09 | 4.20 |  | 5.76E-09 | 4.20 |  | 1.05E-09 | 4.47 |  | 1.05E-09 | 4.49
Table 1: $L_{1}$ norm errors and corresponding convergence rates of the 2D MHD
iso-density vortex (columns 2-5) and the 2D RMHD iso-enthalpy vortex (columns
6-9). Computations were carried out with the resolution of 32, 64, 128, 256
and 512 grid points along each direction.

Both the classical and relativistic cases are evolved until $t=10$, which
corresponds to one period for the MHD case. Table 1 reports the errors for
$B_{x}$ measured in $L_{1}$ norm and the corresponding convergence rates both
for the MHD and RMHD solutions by either using the WENOZ or the MP5
reconstructions.

Fig. 4 overlays the exact pressure profile (continuous black line) on the data
produced with the $4^{\rm th}$-order scheme adopting the WENOZ (red squares),
the MP5 (blue circles) and the $2^{\rm nd}$-order LINEAR scheme (green
triangles). Filled symbols refer to data produced at a resolution of
$N_{x}=128$, while empty ones are produced at $N_{x}=512$ and have been
presented only for the LINEAR scheme, since high-order data do not display a
sensible difference between the resolution of $128$ and $512$. The profiles of
the $4^{\rm th}$-order methods overlap almost perfectly on the exact solution
at $128$, while the $2^{\rm nd}$-order data do not reach such accuracy neither
at $128$ nor at $512$ grid points. Especially in the relativistic case, it is
relevant to notice the significant amount of clipping introduced by the LINEAR
scheme. On the other hand, the $4^{\rm th}$-order method is capable of
attaining significantly improved accuracy already at modest resolutions.

A close comparison with the results of the finite difference (FD) scheme of
Mignone [15] reveals very similar errors for $N_{x}\lesssim 128$ even if the
FD method, being $5^{\rm th}$-order accurate, converges faster than the
proposed FV method. We also point out that, owing to the employment of a
$3^{\rm rd}$-order Runge Kutta method, the time step of the FD scheme has to
be re-scaled as $\Delta t\sim\Delta x^{5/3}$ in order to retain $5^{\rm
th}$-order accuracy (see Eq. [48] of that paper).

###  5.3 Circularly polarized Alfvén waves

We now consider the oblique propagation of circularly polarized Alfvén waves,
which is an exact non-linear solution of the (R)MHD equations and was first
presented by Tóth [6] in the context of classical MHD and later by Del Zanna
[20] for the relativistic case.

In the MHD case we draw upon the setup of [15] by first constructing a 1D wave
profile in the $x$-direction with angular frequency $\omega$, wavenumber $k$
and amplitude $\eta$ as

$\left(\begin{array}[]{c}v_{x}\\\ v_{y}\\\
v_{z}\end{array}\right)=\left(\begin{array}[]{c}0\\\ \eta\sin{\phi}\\\
\eta\cos{\phi}\end{array}\right)\,,\quad\left(\begin{array}[]{c}B_{x}\\\
B_{y}\\\ B_{z}\end{array}\right)=\left(\begin{array}[]{c}c_{a}\sqrt{\rho}\\\
\sqrt{\rho}\eta\sin{\phi}\\\
\sqrt{\rho}\eta\cos{\phi}\end{array}\right)\,,\quad$ (47)

where $\phi=kx-\omega t$, $c_{a}=\omega/k=1$ is the phase (Alfvén) velocity
and $\eta=0.1$ is the wave amplitude. Density and pressure are set to $\rho=1$
and $p=2$, respectively, while we employ the ideal EoS with $\Gamma=5/3$. Such
solution is then rotated by 45∘ around the $y$-axis in the 2D setup and by an
additional 45∘ around the $z$-axis in the 3D case to have the wave propagating
along the computational box’s diagonal. The reader should refer to [15] for
the details.

Figure 5: $L_{1}$-norm errors (left panel) and numerical diffusion (right
panel) for the 3D circularly polarized Alfvén wave test in the MHD regime
(solid lines) and RMHD (dashed-dotted lines) as functions of the grid
resolution $N_{x}$ with $2^{\rm nd}$\- and $4^{\rm th}$-order numerical
schemes. The horizontal dotted line in the left panel is positioned at
$10^{-4}$ in order to highlight the different resolution that achieves an
accuracy of $10^{-4}$ with either a low or a high order scheme. The dashed
lines give the ideal convergence slope (left panel) and a reference slope for
the numerical diffusion (right panel).

The RMHD case is set in a similar way, with a 1D wave profile defined by (see
[20])

$\left(\begin{array}[]{c}v_{x}\\\ v_{y}\\\
v_{z}\end{array}\right)=\left(\begin{array}[]{c}0\\\ -c_{A}\eta\cos{\phi}\\\
-c_{A}\eta\sin{\phi}\end{array}\right)\,,\quad\left(\begin{array}[]{c}B_{x}\\\
B_{y}\\\ B_{z}\end{array}\right)=\left(\begin{array}[]{c}B_{0}\\\
B_{0}\eta\cos{\phi}\\\ B_{0}\eta\sin{\phi}\end{array}\right)\,,\quad$ (48)

where we set $B_{0}=1$, $\eta=1$, $\Gamma=4/3$ and the Alfvén velocity is
defined by

$c_{A}=\frac{B_{0}^{2}}{\rho
h+B_{0}^{2}(1+\eta^{2})}\left[\frac{1}{2}\left(1+\sqrt{1-\left(\frac{2\eta
B_{0}^{2}}{\rho h+B_{0}^{2}(1+\eta^{2})}\right)^{2}}\right)\right]^{-1}\,.$
(49)

We then apply the same procedure as for the MHD case [15] to rotate the wave
front.

The computational domain is the unit square (cube) in 2D (3D), while we follow
the evolution of the wave for one period, after which its profile is supposed
to return exactly on the initial values. We then estimate the error introduced
by the numerical algorithms using the $L_{1}$ norm evaluated for the $y$\-
component of the magnetic field between the final and initial profiles.

Scheme | $N_{x}$ | 2D MP5 |  | 2D WENOZ |  | 3D MP5 |  | 3D WENOZ
---|---|---|---|---|---|---|---|---
|  | $\epsilon_{1}(B_{y})$ | $O_{L_{1}}$ |  | $\epsilon_{1}(B_{y})$ | $O_{L_{1}}$ |  | $\epsilon_{1}(B_{y})$ | $O_{L_{1}}$ |  | $\epsilon_{1}(B_{y})$ | $O_{L_{1}}$
| $8$ | 2.76E-03 | - |  | 3.02E-03 | - |  | 3.61E-03 | - |  | 3.83E-03 | -
| $16$ | 1.47E-04 | 4.23 |  | 1.50E-04 | 4.33 |  | 2.05E-04 | 4.14 |  | 2.09E-04 | 4.19
| $32$ | 8.39E-06 | 4.13 |  | 8.43E-06 | 4.16 |  | 1.24E-05 | 4.05 |  | 1.23E-05 | 4.08
MHD | $64$ | 5.07E-07 | 4.05 |  | 5.08E-07 | 4.05 |  | 7.57E-07 | 4.03 |  | 7.57E-07 | 4.03
| $128$ | 3.13E-08 | 4.02 |  | 3.13E-08 | 4.02 |  | 4.69E-08 | 4.01 |  | 4.69E-08 | 4.01
| $256$ | 1.94E-09 | 4.01 |  | 1.94E-09 | 4.01 |  | 2.92E-09 | 4.01 |  | 2.91E-09 | 4.00
| $512$ | 1.21E-10 | 4.00 |  | 1.21E-10 | 4.00 |  | 1.82E-10 | 4.00 |  | 1.82E-10 | 4.00
| $8$ | 2.27E-02 | - |  | 2.60E-02 | - |  | 3.11E-02 | - |  | 3.48E-02 | -
| $16$ | 1.23E-03 | 4.20 |  | 1.29E-03 | 4.34 |  | 1.75E-03 | 4.15 |  | 1.79E-03 | 4.28
| $32$ | 7.21E-05 | 4.10 |  | 7.25E-05 | 4.15 |  | 1.04E-04 | 4.07 |  | 1.04E-04 | 4.10
RMHD | $64$ | 4.33E-06 | 4.06 |  | 4.33E-06 | 4.05 |  | 6.45E-06 | 4.01 |  | 6.45E-06 | 4.01
| $128$ | 2.67E-07 | 4.02 |  | 2.67E-07 | 4.02 |  | 4.03E-07 | 4.00 |  | 4.03E-07 | 4.00
| $256$ | 1.66E-08 | 4.01 |  | 1.66E-08 | 4.01 |  | 2.52E-08 | 4.00 |  | 2.52E-08 | 3.99
| $512$ | 1.04E-09 | 4.00 |  | 1.04E-09 | 4.00 |  | 1.58E-09 | 4.00 |  | 1.58E-09 | 4.00
Table 2: $L_{1}$ norm errors and corresponding convergence rates computed for
the rotated two-dimensional (columns 2-5) and three-dimensional (columns 6-9)
Alfvén wave propagation. Computations were carried out with the resolution of
8, 16, 32, 64, 128, 256, and 512 grid points along each direction.

Table 2 reports the errors and convergence rates evaluated for a series of 2D
and 3D MHD and RMHD computations employing the MP5 and WENOZ reconstructions
at increasing resolutions $N_{x}=\\{8,16,32,64,128,256,512\\}$. We use a CFL
coefficient of 0.4 and 0.3 for the 2D and 3D case, respectively, and we impose
periodic boundary conditions. The measured convergence rate $O_{L_{1}}$ well
matches the expected value for all the configurations we considered (as shown
in the first panel of Fig. 5). It is interesting to notice that, giving the
same configuration, the $L_{1}$ norm errors are almost identical for the WENOZ
or the MP5 schemes. This confirms that for such a smooth problem, the choice
of the high-order spatial reconstruction algorithm does not change results
significantly.

In the 3D MHD case, our results may be directly compared to those obtained
with the $5^{\rm th}$-order finite difference (FD) scheme of Mignone [15] (see
Table 1 of that paper), indicating that our FV scheme yields smaller errors
for grid sizes lower than $256^{3}$ zones. This difference lessens as the
resolution grows (due to the additional order of accuracy of the FD scheme) so
that similar errors are eventually found at the resolution of $\sim 256^{3}$
zones, where $\epsilon_{1}(B_{y})\sim 3.74\cdot 10^{-9}$ (for the FD scheme)
while $\epsilon_{1}(B_{y})\sim 2.91\cdot 10^{-9}$ (with the present method).

For problems without sharp gradients like the one considered here, the
employment of a $4^{\rm th}$\- (or higher) order scheme leads not only to an
increase in accuracy, but also to a saving in computational time for fixed
accuracy. It is therefore instructive to estimate the net gain of using such
approach. Let $\tau$ be the CPU time for updating a single cell with a chosen
numerical scheme (e.g. $2^{\rm nd}$ or $4^{\rm th}$-order). The total CPU time
$T$ required to update $N^{d}$ cells ($d=1,2,3$ indicates the dimensionality
of problem) scales as

$T=\tau N^{d+1}\,,\quad$ (50)

where the extra factor $N$ comes from the fact that the hyperbolic time step
is inversely proportional to the grid resolution $N$. Now, let
$\epsilon=CN^{-p}$ be the numerical error at the end of the computation, where
$p$ is the order of the scheme and $C$ is an unknown constant depending, in
general, on the chosen numerical method and the problem at hand. A $2^{\rm
nd}$\- and $4^{\rm th}$-order scheme will achieve the same accuracy
$\epsilon_{2}\sim\epsilon_{4}$ using, respectively, $N_{2}$ and $N_{4}$ grid
points related by

$N_{4}\sim\left(\frac{C_{4}}{C_{2}}\right)^{1/4}\sqrt{N_{2}}\,.$ (51)

The previous expression demonstrates that a $4^{\rm th}$-order scheme will
require a reduced grid size of about $\sqrt{N_{2}}$. Furthermore, by
expressing $N$ as a function of $T$ (from Eq. 50) the previous relation can be
adjusted in terms of CPU time:

$\frac{T_{4}}{T_{2}}=\left(\frac{C_{4}}{C_{2}}\right)^{(d+1)/4}\frac{\tau_{4}}{\tau_{2}}\frac{1}{\sqrt{N_{2}^{d+1}}}\,.$
(52)

From our 3D tests we estimated $\tau_{4}/\tau_{2}\sim\\{4.3,3.7\\}$ and
$C_{4}/C_{2}\sim\\{5.0,8.7\\}$ for the MHD and RMHD cases respectively, where
we averaged the results produced with the MP5 and WENOZ spatial
reconstructions. Plugging these parameters into Eq. (52) we obtain a net gain
in computing time, for a given accuracy, that scales as $T_{4}/T_{2}\sim
K/N_{2}^{2}$ with $K\sim\\{21,32\\}$. Our estimates show, for instance, that
for a smooth MHD problem the accuracy obtained with a $2^{\mathrm{nd}}$-order
scheme on a grid with $512\times 512\times 512$ points could be achieved by
our high-order scheme with only $N_{4}\sim 34$ per direction and a computing
time reduced by a factor $\sim 1.2\times 10^{4}$.

Another advantage of using more accurate schemes is the reduction of the
numerical dissipation due to the round-off errors introduced by the grid
discretization. We can estimate the diffusion intrinsic to the numerical
scheme by looking at the damping of Alfvén waves (see [75]). The decay rate
can be expressed as $\Gamma_{A}=\frac{1}{2}\left(\nu+\eta\right)k^{2}$, where
$\nu$ is the fluid’s kinematic viscosity and $\eta$ its resistivity. We then
obtain an approximated value for $\Gamma_{A}$ by computing the quantity
$\Delta\log(\delta B_{z})/\Delta t$, where $\delta B_{z}=\sqrt{\int\int\int
B^{2}_{z}dxdydz}$ and the difference is evaluated after one period. From the
right panel of Fig. 5 we can see how already at very small resolutions ($N\sim
8$) the high-order schemes produce an overall numerical dissipation which is
more than one order of magnitude lower than the one delivered by the
$2^{\mathrm{nd}}$-order scheme. Moreover, the gap between the two schemes
increases significantly at higher resolutions, as the decay rate scales as
$\sim N^{-5}$ in the high-order setup, while we measured a slope of only
$\sim-3$ in the low-order case.

###  5.4 Advection of a magnetic field loop

This problem consists of a weakly magnetized field loop advected by a uniform
velocity field. Since the total pressure is dominated by the thermal
contribution, the magnetic field is essentially transported as a passive
scalar and the preservation of the initial circular shape tests the
dissipative properties of the scheme [9, 71, 14, 15, 76].

Following Gardiner & Stone [9] and Mignone [15], the computational box is
defined by $[-1,1]\times[-0.5,0.5]$ and discretized on $N_{x}\times N_{x}/2$
grid cells. Density and pressure are initially uniformly set to 1 in the MHD
case, while $p=3$ for the RMHD setup (following Beckwith & Stone [76]). The
velocity of the flow is given by

$\mathbf{v}=V_{0}(\cos\alpha,\sin\alpha)\,,$ (53)

with $V_{0}=\sqrt{5}$ in the classical case while $V_{0}=\sqrt{5}/10$ for the
relativistic case. The inclination angle is set to
$\sin\alpha=1/\sqrt{5},\,\cos\alpha=2/\sqrt{5}$. The magnetic field is defined
by means of the vector potential as

$A_{z}=\left\\{\begin{array}[]{ll}A_{0}(R-r)&\textrm{if}\quad R_{1}<r\leq
R\,,\\\ \vskip 6.0pt plus 2.0pt minus 2.0pt\cr 0&\textrm{if}\hskip
9.95863pt\textrm{elsewhere}\,,\end{array}\right.$ (54)

with $A_{0}=10^{-3}$, $R=0.3$, $R_{1}=0.2R$, $r=\sqrt{x^{2}+y^{2}}$.
Simulations have been performed on a computational domain of $256\times 128$
grid points with the LINEAR, WENO3, MP5 and WENOZ schemes up to
$t_{\mathrm{f}}=2$ units for the MHD setup and $t_{\mathrm{f}}=20$ for the
RMHD setup (both corresponding to two advection periods). The CFL was set to
$0.4$ and double periodic boundary conditions are imposed.

Fig. 6 shows the normalized magnetic energy density together with magnetic
field lines with a resolution of $256\times 128$ grid points for each
configuration. As shown in both sequences of Fig. 6, smearing and diffusion
are significantly reduced as we gradually increase the accuracy from $2^{\rm
nd}$ to $4^{\rm th}$-order. It is interesting to notice that the WENOZ
algorithm performs the best in this problem while MP5 has some small amplitude
overshoots (undershoots) at the loop tail (head), most likely caused by a
small clipping of its limiting procedure.

Figure 6: Normalized magnetic energy density for the 2D field loop after two
periods at the resolution of $256\times 128$ grid points in the classical
(top) and relativistic (bottom) MHD. Simulations compare LINEAR (leftmost
panel), WENO3 (center left panels), MP5 (center right panels) and WENOZ
reconstructions (rightmost panels). Magnetic field lines are overlaid as black
lines. Figure 7: Plot of the normalized magnetic energy density decay of a 2D
field loop advected over time at a resolution of $256\times 128$ grid points.
Each simulation has been performed with a $2^{\rm nd}$-order LINEAR (green), a
$3^{\rm rd}$-order WENO3+RK3 (cyan), a $4^{\rm th}$-order MP5 (blue), and a
$4^{\rm th}$-order WENOZ scheme (red). The top panels show the results
obtained in the classical MHD regime, with the HLL (top left panel) and HLLD
(top right panel) Riemann solvers. Similarly, the bottom panels show the
results obtained in the relativistic MHD regime, with the HLL (bottom left
panel) and GFORCE (bottom right panel) Riemann solvers. Notice that for the
results obtained with the HLLD and GFORCE Riemann solvers, the time axes have
been reversed.

Another extremely crucial aspect is the ability of the algorithm to conserve
the energy of the system. Fig. 7 compares the normalized magnetic energy
density over time for a set of $4$ simulations employing a $2^{\rm nd}$\-
(green), a $3^{\rm rd}$\- (cyan), and both the $4^{\rm th}$-order accurate
(blue and red) schemes, respectively, in the MHD (top panels) and RMHD (bottom
panels) regimes. The top panels compare the results of the HLL + UCT-HLL
solver combination with those of the HLLD + UCT-HLLD pair, where the time axis
has been reversed to better compare the two sets of simulations. The RMHD case
follows along but with the GFORCE + UCT-GFORCE solver, showing improved
results, although the choice of the Riemann solver seems to affect only
marginally the conservation of the magnetic energy, with the GFORCE achieving
better performance. The magnetic energy decay is considerably faster for the
$2^{\rm nd}$-order scheme compared to the $3^{\rm rd}$\- and $4^{\rm
th}$-order ones, with the latter yielding the best performance.

Very similar results are also recovered with the FD scheme of Mignone [15]
(see left panel in Fig. A.9 of that paper) when using either WENOZ or MP5,
which dissipate just $\sim 5\%$ of total magnetic energy at the end of the
computation.

###  5.5 2D current sheet

This problem was first introduced by Hawley & Stone [77] and also tested by
Núñez-de La Rosa & Munz both in a MHD [34] and RMHD [35] regime (see also
references therein).

We use the computational box $[-1,1]\times[-0.5,0.5]$ discretized with
$600\times 600$ grid points. Density and pressure of an ideal, adiabatic gas
with $\Gamma=5/3$ are initially set to the constant values $\rho=1$,
$p=0.5\beta$, with $\beta$ being the ratio between the thermal and magnetic
pressure of the gas. The magnetic field initially lies in the $y$-direction
and it is discontinuous with $B_{y}=1$ for $|x|>0.25$ and $B_{y}=-1$
otherwise. We set $v_{x}=\mathrm{Asin(2\pi y)},v_{y}=v_{z}=0$, where $A$ is
the amplitude of the perturbation. The parameters $A$ and $\beta$ are both set
to the challenging values of $A=\beta=0.1$ to test the robustness of the
underlying scheme. Simulations have been run up to $t_{f}=10$, with a CFL
restriction of 0.4 and periodic boundary conditions.

The harmonic velocity perturbation produces non-linear linearly polarized
Alfvén waves that grow into magnetosonic waves, triggering magnetic
reconnection in the regions where the field flips sign. When $\beta<1$,
namely, in the regime of strongly magnetized fluid, magnetic reconnection
forms overpressurized regions that develop magnetosonic waves orthogonal to
the field. This leads to the conversion from magnetic to thermal energy and to
the formation of steep gradients in the regions nearby the field. For this
reason, this benchmark can be extremely severe in assessing the robustness of
the underlying numerical scheme.

In Fig. 8 we display density snapshots at $t=5,7.5,10$ for the MHD (top
panels) and RMHD setups (bottom panels) obtained, respectively, with the HLLC
+ UCT-HLL, and with the GFORCE + UCT-GFORCE as solver plus EMF average
combination. This highly non-linear problem tests our novel limiting strategy
where the WENOZ scheme lowers to WENO3 whenever the order reduction procedure
is activated, thus increasing the stability range of the scheme. A similar
fallback strategy has been adopted by [34, 35] who used a simpler Rusanov-Lax
Friedrichs solver.

Figure 8: Coloured maps for the density of the current sheet test problem for
the classical (top panel) and relativistic (bottom panel) MHD equations. From
left to right snapshots at, respectively, $t=5,7.5,10$ are shown at a
resolution of $600\times 600$ grid points. Simulations are carried out with
the HLLC Riemann solver in the non-relativistic case and with the GFORCE
solver in the relativistic case.

###  5.6 Orszag-Tang vortex

Figure 9: Thermal pressure of the non-relativistic Orszag-Tang vortex computed
at $t=2\pi$ with different algorithms and resolutions. Results are obtained
with a $2^{\rm nd}$-order LINEAR (left column), a $4^{\rm th}$-order MP5
(middle column), and a $4^{\rm th}$-order WENOZ scheme (right column). Each
simulation has been computed with a Roe ($1^{\rm st}$ and $3^{\rm rd}$ rows)
and HLLD ($2^{\rm nd}$ and $4^{\rm th}$ rows) Riemann solver with a grid
resolution of $N_{x}=N_{y}=64$ ($1^{\rm st}$ and $2^{\rm nd}$ rows) and
$N_{x}=N_{y}=512$ ($3^{\rm rd}$ and $4^{\rm th}$ rows). Figure 10: Thermal
pressure of the relativistic Orszag-Tang vortex computed at $t=2\pi$ with
different algorithms and resolutions. From left to right, results are obtained
with, respectively, a $2^{\rm nd}$-order LINEAR (left column), a $4^{\rm
th}$-order MP5 (middle column), and a $4^{\rm th}$-order WENOZ scheme (bottom
column). The grid resolution adopted is $N_{x}=N_{y}=128$ (top row) and
$N_{x}=N_{y}=1024$ (bottom row). The presence or absence of O-points emerging
from the central current sheets (only for the high resolution cases) is
highlighted by a zoom on the central region.

The Orszag-Tang MHD vortex [78] has traditionally been a standard numerical
benchmark for inter-scheme comparison. The test consists of an initial
vorticity distribution that spins the fluid counter-clockwise, leading
eventually to the formation of a complex structure that includes shocks and
turbulence. The initial density and pressure are set to the constant values
$\rho_{0}=25/9$ and $p_{0}=5/3$, while the initial velocity and magnetic
fields are $(v_{x},v_{y})=(-\sin{y},\sin{x})$ and
$(B_{x},B_{y})=(-\sin{y},\sin{2x})$, respectively, in the MHD setup, and
$(v_{x},v_{y})=(-0.5\sin{y},0.5\sin{x})$ in the RMHD setup (as presented in
[79, 80, 76]). We ran computations up to $t_{f}=2\pi$ with
$N_{x}=N_{y}=\\{64,512\\}$ grid zones for the MHD case, and
$N_{x}=N_{y}=\\{128,1024\\}$ for the RMHD case using periodic boundary
conditions and a CFL number of $0.5$. We repeated the computations by choosing
two different Riemann solvers in the MHD case (Roe and HLLD, always with the
UCT-HLLD average), and only the GFORCE Riemann solver coupled to the UCT-
GFORCE EMF averaging algorithm in the relativistic case. Likewise, we also
compared 3 different reconstruction schemes, namely, LINEAR, MP5, and WENOZ.

Multiple shock-vortex interactions regulate the dynamics of the system,
leading to the formation of an inclined current sheet at the center of the
domain. The magnetic energy is gradually dissipated and the current sheet
twists, leading to the structures observed in the thermal pressure
distribution shown in Fig. 9. For the MHD setup, even at modest resolutions
(64 points) the $4^{\rm th}$-order scheme with MP5 reconstruction and the HLLD
Riemann solver allows to distinguish at the center of the computational domain
the formation of a magnetic island (an O-point). It is well known (see [81,
82]) that the tearing instability is expected to develop on the ideal
(Alfvénic) timescales of an ideal MHD simulation only if the numerical
dissipation is sufficiently low. At higher resolutions, all simulations
produce high-pressure regions, albeit small-scale structures are best
appreciated with the high-order runs.

Similarly to Mignone & Del Zanna [18], we can assess more quantitatively the
improvements introduced by our $4^{\rm th}$-order scheme by computing the
parameter $\mu_{p}=\mathrm{max}(p)/\langle p\rangle$, i.e. the pressure
maximum normalized by the average pressure at the end of a run. As this
quantity generally increases with better resolution, we use it as a proxy to
assess the decrease of numerical dissipation. For instance, the values of
$\mu_{p}$ obtained with the MP5 reconstruction and HLLD solver against those
obtained with the LINEAR one yield an average increase in effective resolution
of $\sim 2.23$, which leads to an average speedup of $\sim 2.58$ (taking into
account the extra computational cost of the $4^{\rm th}$-order scheme
estimated in §5.3). This is, however, only a rough estimate that cannot fully
characterize the impact of our scheme on the typical numerical dissipation
length-scale in a turbulent plasma, for which more dedicated tests and
analysis are required.

For the RMHD case, the final stage of the simulation produces a central low-
pressure cavity with a thin high-pressure filament, which exhibits finer
structures for less dissipative schemes. While at a resolution of $128$ points
there are no apparent deviations between the different setups, at higher
resolution ($1024$ points) the numerical resistivity from the $4^{\rm
th}$-order scheme is low enough to lead to the production of multiple O-points
in the high-pressure filament, although they remain still quite unresolved.
Our results confirm the overall robustness of the scheme as well as its
ability to capture small-scale structures that would otherwise be lost with a
lower order scheme.

###  5.7 3D Blast Wave

Case | $\rho_{\rm in}$ | $\rho_{\rm out}$ | $p_{\rm in}$ | $p_{\rm out}$ | $B_{0}$ | $\Gamma$ | $r_{0}$ | $L$ | $t_{\rm stop}$ | Solver | En. Corr.
---|---|---|---|---|---|---|---|---|---|---|---
C1 (MHD) | $1$ | $1$ | $10$ | $0.1$ | $3$ | $5/3$ | $0.1$ | $1$ | $0.2$ | HLL | NO
C2 (MHD) | $1$ | $1$ | $10^{3}$ | $0.1$ | $100/\sqrt{4\pi}$ | $1.4$ | $0.1$ | $1$ | $0.01$ | HLL | YES
R1 (RMHD) | $10^{-2}$ | $10^{-4}$ | $1$ | $5\times 10^{-3}$ | $1$ | $4/3$ | $0.08$ | $6$ | $4$ | GFORCE | YES
Table 3: Parameters used in the 3D blast wave problem.

We finally consider the 3D magnetized blast wave problem in order to assess
the robustness of the algorithms under strong magnetization conditions.
Despite its simplicity, the blast wave problem is a particularly effective
benchmark in testing the solver’s ability at handling MHD wave degeneracies
parallel and perpendicularly to the field orientation.

The initial condition consists of a uniform medium with density and pressure
values set, respectively, to $\rho_{\rm out}$ and $p_{\rm out}$ and threaded
by an oblique constant magnetic field,

$\mathbf{B}=B_{0}\Big{(}\sin{\theta}\cos{\phi},\,\sin{\theta}\sin{\phi},\,\cos{\theta}\Big{)}$
(55)

where $\theta=\pi/2$ and $\phi=\pi/4$ are used in what follows. Note also that
oblique configurations make the computation considerably more challenging than
in the grid-aligned cases (see, for instance the discussion in §4.4 of Mattia
& Mignone [83]). The computational domain is a cube of length $L$ centered at
the origin with zero-gradient boundary conditions. An over-pressurized
spherical region where pressure and density can have different values, namely
$\rho_{\rm in}$ and $p_{\rm in}$, is initially set for $r<r_{0}$.

We examine three different cases with parameters listed in table 3. Cases C1
is taken from Felker & Stone [36], with the exception that we consider a
magnetization three times as large here. Case C2 proposes a more severe
configuration as initially presented by Balsara & Spicer [5] (see also [18]).
Lastly, in case R1 (from [83]) we apply our $4^{\rm th}$-order method to the
solution of the relativistic MHD equations. Note that the energy correction
step (see [5] and, e.g., see Sec. 6.4 of [18]) is applied to both cases C2 and
R1, while case C1 has the largest magnetization our scheme can afford without
it. This correction step is still performed at $2^{\rm nd}$-order level and it
was found to substantially improve the robustness of the code for larger
magnetic field strengths. Computations are performed using $192\times
192\times 192$ grid zones.

Figure 11: MHD 3D blast wave problem (case C1) at $t=0.2$. The left and
central panels show density (top left), pressure (top middle, in logarithmic
scale), specific kinetic energy (bottom left) and magnetic energy (bottom
right) at $z=0$ for a simulation employing $192\times 192\times 192$ grid
points. In the rightmost panels we show the density (top) and magnetic energy
density (bottom) on the major (red line) and minor (black line) diagonals.
Figure 12: Same as Fig. 11 but for the case C2 at $t=0.01$. Figure 13: Same as
Fig. 11 but for the case R1 at $t=4$. In the bottom left panel, the specific
kinetic energy has been replaced by the relativistic counterpart.

Results are shown in Fig. 11 (C1), 12 (C2) and 13 (R1) for the three selected
cases. The explosion spawns an oval-shaped structure delimited by an outer
fast forward shock and the presence of a magnetic field makes the propagation
highly anisotropic by compressing the gas in the direction parallel to the
field. Gas motion takes place along the field lines generating two oppositely
moving plasma blobs. In the perpendicular direction, the outer fast shock
becomes magnetically dominated with very weak compression. Discontinuous waves
and smooth structures are well resolved by our scheme and the point-symmetry
is preserved throughout the evolution. We point out that the fallback approach
played an essential role for this test thus strengthening the benefits of a
low/high order hybrid scheme while resolving complex flows featuring both
smooth as well as discontinuous solutions.

###  5.8 3D Cloud-Shock Interaction

We propose a 3D version of the cloud-shock interaction test previously
introduced by [72, 84, 85, 29] in its MHD version, and by Mignone [69] for the
RMHD case. For both configurations, the initial condition consists of a plane
discontinuity propagating a fast shock and a rotational wave. Specifically, in
the MHD case, the initial state
$\altmathcal{V}=(\rho,v_{x},v_{y},v_{z},B_{x},B_{y},B_{z},p)$ is initialized
as

$\left\\{\begin{array}[]{lclr}\altmathcal{V}_{L}&=&\displaystyle\left(3.86859,\,0,\,0,\,0,\,0,\,2.1826182,\,-2.1826182,\,167.345\right)&\quad{\rm
for}\quad x<0.6\,,\\\ \vskip 6.0pt plus 2.0pt minus
2.0pt\cr\altmathcal{V}_{R}&=&\displaystyle\left(1,\,-11.2536,\,0,\,0,\,0,\,0.56418958,\,0.56418958,\,1\right)&{\rm
for}\quad x>0.6\,,\end{array}\right.$ (56)

while in the RMHD test as

$\left\\{\begin{array}[]{lclr}\altmathcal{V}_{L}&=&\displaystyle\left(39.5052,\,0,\,0,\,0,\,0,\,0,\,1.9753,\,129.72386\right)&\quad{\rm
for}\quad x<0.6\,,\\\ \vskip 6.0pt plus 2.0pt minus
2.0pt\cr\altmathcal{V}_{R}&=&\displaystyle(1,\,-\sqrt{0.99},\,0,\,0,\,0,\,0,\,0.5,\,10^{-3})&{\rm
for}\quad x>0.6\,.\end{array}\right.$ (57)

At $\vec{x}=(0.8,0.5,0.5)$ an overdense ($\rho=10$) spherical cloud of radius
$0.15$ is initially in hydrostatic equilibrium with the ambient plasma under
the adiabatic equation of state with $\Gamma=5/3$ ($\Gamma=4/3$ in the special
relativistic case). For both setups the computational domain consists of the
unit cube, with outflow (zero-gradient) boundary conditions employed
everywhere, except at the lower $x$\- boundary, where the initial condition is
imposed.

During the evolution, the magnetic shock propagates to the right, and the high
density cloud moves supersonically into the shock front. At the final
evolution time ($t_{f}=0.06$ MHD, $t_{f}=0.6$ RMHD) the cloud has been
completely shocked and it has developed its characteristic mushroom-like
structure. Fig. 14 shows the initial and final density condition for a
$256\times 256\times 256$ grid points simulation. The $4^{\rm th}$-order
scheme properly recovers the expected structures. It is important to mention
that, for this particular test problem, the employment of the order reduction
procedure extended to the selective fluxes integration has been fundamental
for the correct evolution of the simulations.

Figure 14: Non-relativistic (top panels) and relativistic (bottom panels)
cloud-shock interaction at initial ($t=0$, left panels) and final time
($t=0.06$ for the non-relativistic case, and $t=0.6$ for the relativistic
case, right panels). Shown is the density at $z=0.5$ in logarithmic scale.
Both simulations ran with $256\times 256\times 256$ grid points. The
colorscales span over a range of $[1,24]$ for the non-relativistic case and
$[1,110]$ for the relativistic case.

##  6 Conclusions

In this paper we have presented a genuinely $4^{\rm th}$-order accurate finite
volume method for the solution of the classical and special relativistic MHD
equations employing the constrained transport algorithm. Our scheme has been
implemented in the PLUTO code for astrophysical plasma dynamics [68] and is
rooted over the method originally proposed by McCorquodale & Colella [28],
adding a number of innovative aspects that yield an accurate, robust and
efficient computational tool. In fact, our scheme maximizes the benefits of
the FV formalism (e.g., robustness when dealing with shocks, discontinuities,
and steep gradients) along with techniques inspired by the FD formalism, such
as the introduction of pointwise to pointwise reconstruction operations that
ease up the structure of the scheme. To this aim, we revisited standard
spatial reconstruction schemes (such as WENOZ and MP5), traditionally based on
the knowledge of 1D volume averages, so as to provide $4^{\rm th}$-order
accurate left and right states at zone interfaces by directly interpolating
point values rather than volume averages. Pointwise reconstructions can be
accomplished in primitive or characteristic variables and they lead to a
compact stencil and a more cost-effective scheme. Moreover, our implementation
is endowed with the generalized upwind constrained transport (UCT) of Mignone
& Del Zanna [18], plus a new generalization of the UCT-GFORCE averaging scheme
for special relativistic MHD that has been widely validated in the presented
benchmarks.

The performance of the method has been assessed against selected 2D and 3D
numerical test problems, demonstrating: i) genuine $4^{\rm th}$-order accuracy
in smooth problems, ii) reduced numerical dissipation iii) the code’s ability
to correctly capture propagating discontinuities and to resolve small-scale
structures and turbulence. The robustness of the method is reinforced by a
fallback approach that reverts the spatial order of the scheme to $2$ in
proximity of steep gradients. To identify such critical regions, we have
proposed a new discontinuity detector - based on a derivatives ratio - which
employs up to the fourth derivative to pinpoint high-frequency oscillations
typically arising in proximity of a discontinuity. The order reduction
procedure turned out to be, in our experience, a crucial element when
undertaking severe configurations, allowing our high-order scheme to
successfully carry out computations otherwise feasible only with $2^{\rm
nd}$-order schemes (see, e.g., configuration C2 §5.7).

We point out that, despite the augmented complexity of high-order algorithms,
genuine $4^{\rm th}$-order accuracy in smooth problems leads to a net saving
of computational time. In §5.3 we demonstrated that, for a given accuracy, the
high-order scheme will require, approximately $\sqrt{N_{2}}$ grid cells when
compared to a $2^{\rm nd}$-order scheme with a corresponding reduction of CPU
time that scales as $\sim N_{2}^{-(d+1)/2}$, where $d$ is the number of
spatial dimensions. When applied to simple 3D computations, for instance, we
found that the accuracy reached by a $2^{\mathrm{nd}}$-order scheme with
$512^{3}$ grid points is equally obtained when employing our $4^{\rm
th}$-order scheme with only $\sim 34^{3}$ points, with a corresponding
reduction of the total computational time by a few orders of magnitude.

Better accuracy leads also to reduced numerical dissipation and thus an
enhanced potential at describing spatial scales that would otherwise be
challenging to probe within a traditional $2^{\rm nd}$-order framework. Our
results indicate that the decay rate of numerical diffusion is $\sim N^{-5}$
for our novel scheme, and only $\sim N^{-3}$ for the traditional scheme. For
all these reasons our $4^{\rm th}$-order scheme proved to be a computational
tool suited for the next generation of numerical simulations in computational
fluid dynamics. In a companion paper we will extend this formalism to non-
Cartesian geometries.

Compared to other high-order competitive methods, such as discontinuous
Galerkin (DG), our FV scheme requires fewer equations to be evolved in time
(e.g. $8$ instead of $32$ for a $4^{\rm th}$-order method), does not reduce
stability with increasing order and it is better suited for flows containing
smooth and discontinuous features. Its parallel implementation requires,
however, several more inter-processor communications (one per Runge-Kutta
stage). A more thorough comparison, though, has to be considered in a separate
paper.

Acknowledgments. We thank L. Del Zanna for useful discussions and comments.
The PLUTO code is publicly available and the simulation data will be shared on
reasonable request to the corresponding author.

This project has received funding from the European Union’s Horizon Europe
research and innovation programme under the Marie Skłodowska-Curie grant
agreement No 101064953 (GR-PLUTO).

This work has received funding from the European High Performance Computing
Joint Undertaking (JU) and Belgium, Czech Republic, France, Germany, Greece,
Italy, Norway, and Spain under grant agreement No 101093441 (SPACE).

This paper is supported by the Fondazione ICSC, Spoke 3 Astrophysics and
Cosmos Observations. National Recovery and Resilience Plan (Piano Nazionale di
Ripresa e Resilienza, PNRR) Project ID CN_00000013 “Italian Research Center on
High-Performance Computing, Big Data and Quantum Computing” funded by MUR
Missione 4 Componente 2 Investimento 1.4: Potenziamento strutture di ricerca e
creazione di “campioni nazionali di $R\&S$ (M4C2-19 )” - Next Generation EU
(NGEU).

We also acknowledge the OCCAM supercomputing facility available at the
Competence Centre for Scientific Computing at the University of Torino.

##  Appendix A Laplacian Operators in Cartesian Coordinates

We start by demonstrating the validity of Eq. (16) and deriving the expression
of the $2^{\rm nd}$-order accurate Laplacian operator for a given function
$U(x)$ in a Cartesian mesh in one dimension. The possibility to
straightforwardly generalize the treatment to multiple dimensions is
guaranteed when considering Cartesian meshes due to the orthogonality of the
standard basis. The (1D) cell-average is given by

$\displaystyle\left<{U}\right>_{{\boldsymbol{c}}}=\frac{1}{h}\int_{-\frac{h}{2}}^{\frac{h}{2}}U(x)\,dx\,,$
(58)

where $h$ is the cell width. By expanding in Taylor series $U(x)$ near the
cell-center $x_{c}$ up to $4^{\rm th}$-order we obtain

$\displaystyle\left<{U}\right>_{{\boldsymbol{c}}}=\frac{1}{h}\int_{-\frac{h}{2}}^{\frac{h}{2}}\left[U(x_{c})+U^{\prime}(x_{c})(x-x_{c})+\frac{U^{\prime\prime}(x_{c})}{2}(x-x_{c})^{2}+\frac{U^{\prime\prime\prime}(x_{c})}{6}(x-x_{c})^{3}+O(x^{4})\right]\,dx\,,$
(59)

and, since odd integrands vanish when integrated over symmetrical domains,
from straightforward calculation we get

$\displaystyle\left<{U}\right>_{{\boldsymbol{c}}}=U(x_{c})+\frac{h^{2}}{24}U^{\prime\prime}(x_{c})+O(h^{4})\,.$
(60)

We now obtain a $2^{\rm nd}$-order accurate expression of the $2^{\rm
nd}$-order derivative $U^{\prime\prime}$. By expanding in Taylor series up to
$4^{\rm th}$-order the nearby average values:
$\left<{U}\right>_{{\boldsymbol{c}}-\hat{\mathbf{e}}_{x}}$,
$\left<{U}\right>_{{\boldsymbol{c}}+\hat{\mathbf{e}}_{x}}$ we get

$\begin{array}[]{lc}\displaystyle\left<{U}\right>_{{\boldsymbol{c}}-\hat{\mathbf{e}}_{x}}=\frac{1}{h}\int_{-\frac{3h}{2}}^{-\frac{h}{2}}\left[U(x_{c})+U^{\prime}(x_{c})(x-x_{c})+\frac{U^{\prime\prime}(x_{c})}{2}(x-x_{c})^{2}+\frac{U^{\prime\prime\prime}(x_{c})}{6}(x-x_{c})^{3}+O(x^{4})\right]\,dx\,,\\\
\vskip 6.0pt plus 2.0pt minus
2.0pt\cr\displaystyle\left<{U}\right>_{{\boldsymbol{c}}+\hat{\mathbf{e}}_{x}}=\frac{1}{h}\int_{\frac{h}{2}}^{\frac{3h}{2}}\left[U(x_{c})+U^{\prime}(x_{c})(x-x_{c})+\frac{U^{\prime\prime}(x_{c})}{2}(x-x_{c})^{2}+\frac{U^{\prime\prime\prime}(x_{c})}{6}(x-x_{c})^{3}+O(x^{4})\right]\,dx\,.\end{array}$
(61)

Since the interval is no longer symmetric, the odd powers no longer integrate
to zero. The final result is

$\displaystyle\left<{U}\right>_{{\boldsymbol{c}}-\hat{\mathbf{e}}_{x}}$
$\displaystyle=U(x_{c})-U^{\prime}(x_{c})h+\frac{13}{24}U^{\prime\prime}(x_{c})h^{2}-\frac{5}{24}U^{\prime\prime\prime}(x_{c})h^{3}+O(h^{4})\,,$
(62) $\displaystyle\left<{U}\right>_{{\boldsymbol{c}}+\hat{\mathbf{e}}_{x}}$
$\displaystyle=U(x_{c})+U^{\prime}(x_{c})h+\frac{13}{24}U^{\prime\prime}(x_{c})h^{2}+\frac{5}{24}U^{\prime\prime\prime}(x_{c})h^{3}+O(h^{4})\,.$
(63)

Adding Eq. (62) to Eq. (63) and subtracting twice Eq. (60)

$\displaystyle\left<{U}\right>_{{\boldsymbol{c}}+\hat{\mathbf{e}}_{x}}+\left<{U}\right>_{{\boldsymbol{c}}-\hat{\mathbf{e}}_{x}}-2\left<{U}\right>_{{\boldsymbol{c}}}=U^{\prime\prime}h^{2}+O(h^{4})\,,$
(64)

we can express the $2^{\rm nd}$-order accurate second derivative in terms of
volume-averaged quantities as

$\displaystyle
U^{\prime\prime}=\frac{\left<{U}\right>_{{\boldsymbol{c}}+\hat{\mathbf{e}}_{x}}-2\left<{U}\right>_{{\boldsymbol{c}}}+\left<{U}\right>_{{\boldsymbol{c}}-\hat{\mathbf{e}}_{x}}}{h^{2}}+O(h^{2})\,,$
(65)

which defines the Laplacian operator in Eq. (16).

Conversely, Eq. (17) can be formally obtained applying a Simpson’s quadrature
rule by integrating over the cell width $h$

$\displaystyle\int_{-\frac{h}{2}}^{\frac{h}{2}}U(x)\,dx\simeq\displaystyle\int_{-\frac{h}{2}}^{\frac{h}{2}}\left(\alpha
x^{2}+\beta x+\gamma\right)\,dx\,.$ (66)

The function $U(x)$ is approximated by a $2^{\rm nd}$-order polynomial in the
coefficients $\alpha,\beta,\gamma$. The polynomial in the right hand side of
Eq. (66) is the best approximation of $U(x)$ when the coefficients
$\alpha,\beta,\gamma$ are set by evaluating the function in the $3$ nearby
cell centers:

$\displaystyle\alpha=\displaystyle\frac{U_{{\boldsymbol{c}}+\hat{\mathbf{e}}_{x}}-2U_{{\boldsymbol{c}}}+U_{{\boldsymbol{c}}-\hat{\mathbf{e}}_{x}}}{2}\,,\quad\displaystyle\beta=\displaystyle\frac{U_{{\boldsymbol{c}}+\hat{\mathbf{e}}_{x}}-U_{{\boldsymbol{c}}-\hat{\mathbf{e}}_{x}}}{2}\,,\quad\displaystyle\gamma=\displaystyle
U_{{\boldsymbol{c}}}\,.$ (67)

By replacing these expressions in Eq. (66) and by direct calculation of the
integral we obtain Eq. (17).

##  Appendix B Explicit Expressions for the UCT Scheme

Here, for the sake of completeness and to also amend some equations containing
typos, we briefly recap the most relevant expressions needed in our UCT
scheme, as originally presented by Mignone & Del Zanna [18]. The edge-centered
EMF (Eq. 39) requires the reconstructed values of the transverse velocities as
well as the linear combination coefficients $a$ and $d$ (for the centred and
diffusive flux contribution, respectively).

  * •

Transverse Velocities: when using the HLL, HLLC or HLLD classical MHD Riemann
solvers at a zone interface, the transverse velocities are obtained, e.g., at
an $x$-face, as

$\overline{v}_{t,{\mathbf{x}_{f}}}=\frac{\alpha_{x}^{R}v^{L}_{t,{\mathbf{x}_{f}}}+\alpha_{x}^{L}v^{R}_{t,{\mathbf{x}_{f}}}}{\alpha^{R}_{x}+\alpha^{L}_{x}}\,,\quad(t=y,z)\,,$
(68)

where $v^{L/R}_{t,{\mathbf{x}_{f}}}$ are the left / right transverse
velocities reconstructed at an $x$-face,

$\alpha_{x}^{R}=\max(0,\lambda^{R}_{{\mathbf{x}_{f}}}),\quad\alpha_{x}^{L}=-\min(0,\lambda^{L}_{{\mathbf{x}_{f}}})$
(69)

with $\lambda^{L/R}_{{\mathbf{x}_{f}}}$ being the smallest and largest
characteristic speeds (see Eq. 14 of [18]). In the case of the GFORCE Riemann
solver instead, we use

$\overline{v}_{t}=\omega\left[v_{t}^{\rm
LW}-\frac{\tau_{x}}{2}\left(v^{R}_{t}-v_{t}^{L}\right)v_{x}^{\rm
LW}\right]+(1-\omega)\left(\frac{v_{t}^{L}+v_{t}^{R}}{2}\right)$ (70)

(note that Eq. 51 of [18] contains “$y$” instead of “$t$”), with
$\omega=1/(1+c)$ ($c$ is the Courant factor, $c=1$ yields the FORCE scheme),
$\mathbf{v}^{\rm LW}$ is the velocity in the Lax-Wendroff state (see. Eq. 48
of [18]), and
$\tau_{x}=1/\max(|\lambda^{L}_{{\mathbf{x}_{f}}}|,\,|\lambda^{R}_{{\mathbf{x}_{f}}}|)$.

  * •

Combination Coefficients $a$ and $d$. For the HLL or HLLC Riemann solvers
(away from degeneracies), the combination coefficients are evaluated as
arithmetic averages across two adjacent faces (Eq. 34 and 35 of [18]) using
the values available with the 1D Riemann solver:

$a^{L}_{x}=\frac{\alpha_{x}^{R}}{\alpha^{R}_{x}+\alpha^{L}_{x}}\,,\quad
a^{R}_{x}=\frac{\alpha_{x}^{L}}{\alpha^{R}_{x}+\alpha^{L}_{x}}\,,\quad
d^{L}_{x}=d^{R}_{x}=\frac{\alpha^{R}_{x}\alpha^{L}_{x}}{\alpha^{R}_{x}+\alpha^{L}_{x}}\,,$
(71)

where $\alpha_{x}^{L/R}$ have been given in Eq. (69). In the limit $B_{x}\to
0$, the corresponding expressions for the HLLC solver are given by Eq. 38 of
[18]. For the HLLD solver, instead, we employ

$a^{L}=\frac{1+\nu^{*}}{2}\,,\quad a^{R}=\frac{1-\nu^{*}}{2}\,,\quad
d^{s}=\frac{1}{2}(\nu^{s}-\nu^{*})\tilde{\chi}^{s}+\frac{1}{2}\left(|\lambda^{*s}|-\nu^{*}\lambda^{*s}\right)\,,\quad({\rm
s}=L,R)$ (72)

where $\tilde{\chi}^{s}=(\lambda^{*s}-\lambda^{s})\chi^{s}$, while

$\nu^{s}=\frac{\lambda^{*s}+\lambda^{s}}{|\lambda^{*s}|+|\lambda^{s}|}\,,\qquad\nu^{*}=\frac{\lambda^{*R}+\lambda^{*L}}{|\lambda^{*R}|+|\lambda^{*L}|}\,,$
(73)

and $\lambda^{*L/R}$ are the Alfvén velocities of the HLLD fan (see Eq. 40 of
[18]).

## References

  * [1] S. A. E. G. Falle, S. S. Komissarov, P. Joarder, A multidimensional upwind scheme for magnetohydrodynamics, MNRAS297 (1) (1998) 265–277.  doi:10.1046/j.1365-8711.1998.01506.x.
  * [2] D. S. Balsara, Total Variation Diminishing Scheme for Adiabatic and Isothermal Magnetohydrodynamics, ApJS116 (1) (1998) 133–153.  doi:10.1086/313093.
  * [3] W. Dai, P. R. Woodward, On the Divergence-free Condition and Conservation Laws in Numerical Simulations for Supersonic Magnetohydrodynamical Flows, ApJ494 (1) (1998) 317–335.  doi:10.1086/305176.
  * [4] K. G. Powell, P. L. Roe, T. J. Linde, T. I. Gombosi, D. L. De Zeeuw, A Solution-Adaptive Upwind Scheme for Ideal Magnetohydrodynamics, Journal of Computational Physics 154 (2) (1999) 284–309.  doi:10.1006/jcph.1999.6299.
  * [5] D. S. Balsara, D. S. Spicer, A Staggered Mesh Algorithm Using High Order Godunov Fluxes to Ensure Solenoidal Magnetic Fields in Magnetohydrodynamic Simulations, Journal of Computational Physics 149 (1999) 270–292.  doi:10.1006/jcph.1998.6153.
  * [6] G. Tóth, The $\nabla$· B=0 Constraint in Shock-Capturing Magnetohydrodynamics Codes, Journal of Computational Physics 161 (2) (2000) 605–652.  doi:10.1006/jcph.2000.6519.
  * [7] D. S. Balsara, Second-Order-accurate Schemes for Magnetohydrodynamics with Divergence-free Reconstruction, ApJS151 (1) (2004) 149–184.  arXiv:astro-ph/0308249, doi:10.1086/381377.
  * [8] M. Torrilhon, Locally divergence-preserving upwind finite volume schemes for magnetohydrodynamic equations, SIAM J. Scientific Computing 26 (2005) 1166–1191.  doi:10.1137/S1064827503426401.
  * [9] T. A. Gardiner, J. M. Stone, An unsplit Godunov method for ideal MHD via constrained transport, Journal of Computational Physics 205 (2005) 509–539.  arXiv:arXiv:astro-ph/0501557, doi:10.1016/j.jcp.2004.11.016.
  * [10] S. Fromang, P. Hennebelle, R. Teyssier, A high order Godunov scheme with constrained transport and adaptive mesh refinement for astrophysical magnetohydrodynamics, A&A457 (2) (2006) 371–384.  arXiv:astro-ph/0607230, doi:10.1051/0004-6361:20065371.
  * [11] A. Mignone, G. Bodo, An HLLC Riemann solver for relativistic flows - II. Magnetohydrodynamics, MNRAS368 (35) (2006) 1040–1054.  arXiv:astro-ph/0601640, doi:10.1111/j.1365-2966.2006.10162.x.
  * [12] J. Rossmanith, An unstaggered, high?resolution constrained transport method for magnetohydrodynamic flows, SIAM J. Scientific Computing 28 (2006) 1766–1797.  doi:10.1137/050627022.
  * [13] G. Tóth, D. L. De Zeeuw, T. I. Gombosi, K. G. Powell, A parallel explicit/implicit time stepping scheme on block-adaptive grids, Journal of Computational Physics 217 (2) (2006) 722–758.  doi:10.1016/j.jcp.2006.01.029.
  * [14] D. Lee, A. E. Deane, An unsplit staggered mesh scheme for multidimensional magnetohydrodynamics, Journal of Computational Physics 228 (4) (2009) 952–975.  doi:10.1016/j.jcp.2008.08.026.
  * [15] A. Mignone, P. Tzeferacos, G. Bodo, High-order conservative finite difference GLM-MHD schemes for cell-centered MHD, Journal of Computational Physics 229 (17) (2010) 5896–5920.  arXiv:1001.2832, doi:10.1016/j.jcp.2010.04.013.
  * [16] C. D. Meyer, D. S. Balsara, T. D. Aslam, A second-order accurate Super TimeStepping formulation for anisotropic thermal conduction, MNRAS422 (3) (2012) 2102–2115.  doi:10.1111/j.1365-2966.2012.20744.x.
  * [17] A. Mignone, M. Flock, M. Stute, S. M. Kolb, G. Muscianisi, A conservative orbital advection scheme for simulations of magnetized shear flows with the PLUTO code, A&A545 (2012) A152.  arXiv:1207.2955, doi:10.1051/0004-6361/201219557.
  * [18] A. Mignone, L. Del Zanna, Systematic construction of upwind constrained transport schemes for MHD, Journal of Computational Physics 424 (2021) 109748\.  arXiv:2004.10542, doi:10.1016/j.jcp.2020.109748.
  * [19] P. Londrillo, L. Del Zanna, High-Order Upwind Schemes for Multidimensional Magnetohydrodynamics, ApJ530 (1) (2000) 508–524.  arXiv:astro-ph/9910086, doi:10.1086/308344.
  * [20] L. Del Zanna, O. Zanotti, N. Bucciantini, P. Londrillo, ECHO: a Eulerian conservative high-order scheme for general relativistic magnetohydrodynamics and magnetodynamics, A&A473 (2007) 11–30.  arXiv:0704.3206, doi:10.1051/0004-6361:20077093.
  * [21] A. J. Christlieb, J. A. Rossmanith, Q. Tang, Finite difference weighted essentially non-oscillatory schemes with constrained transport for ideal magnetohydrodynamics, Journal of Computational Physics 268 (2014) 302–325.  arXiv:1309.3344, doi:10.1016/j.jcp.2014.03.001.
  * [22] Y. Chen, G. Tóth, T. I. Gombosi, A fifth-order finite difference scheme for hyperbolic equations on block-adaptive curvilinear grids, Journal of Computational Physics 305 (2016) 604–621.  doi:10.1016/j.jcp.2015.11.003.
  * [23] T. Minoshima, T. Miyoshi, Y. Matsumoto, A High-order Weighted Finite Difference Scheme with a Multistate Approximate Riemann Solver for Divergence-free Magnetohydrodynamic Simulations, ApJS242 (2) (2019) 14.  arXiv:1903.04759, doi:10.3847/1538-4365/ab1a36.
  * [24] A. Reyes, D. Lee, C. Graziani, P. Tzeferacos, A variable high-order shock-capturing finite difference method with GP-WENO, Journal of Computational Physics 381 (2019) 189–217.  arXiv:1808.04960, doi:10.1016/j.jcp.2018.12.028.
  * [25] C.-W. Shu, High Order Weighted Essentially Nonoscillatory Schemes for Convection Dominated Problems, SIAM Review 51 (1) (2009) 82–126.  doi:10.1137/070679065.
  * [26] D. S. Balsara, Divergence-free reconstruction of magnetic fields and WENO schemes for magnetohydrodynamics, Journal of Computational Physics 228 (14) (2009) 5040–5056.  arXiv:0811.2192, doi:10.1016/j.jcp.2009.03.038.
  * [27] P. Colella, M. R. Dorr, J. A. F. Hittinger, D. F. Martin, High-order, finite-volume methods in mapped coordinates, Journal of Computational Physics 230 (8) (2011) 2952–2976.  doi:10.1016/j.jcp.2010.12.044.
  * [28] P. McCorquodale, P. Colella, A high-order finite-volume method for conservation laws on locally refined grids., Commun. Appl. Math. Comput. Sci. 6 (1) (2011) 1–25.  doi:10.2140/camcos.2011.6.1.
  * [29] C. Helzel, J. A. Rossmanith, B. Taetz, An unstaggered constrained transport method for the 3D ideal magnetohydrodynamic equations, Journal of Computational Physics 230 (10) (2011) 3803–3829.  arXiv:1007.2606, doi:10.1016/j.jcp.2011.02.009.
  * [30] A. Susanto, L. Ivan, H. De Sterck, C. P. T. Groth, High-order central ENO finite-volume scheme for ideal MHD, Journal of Computational Physics 250 (2013) 141–164.  doi:10.1016/j.jcp.2013.04.040.
  * [31] D. S. Balsara, Higher-order accurate space-time schemes for computational astrophysics—Part I: finite volume methods, Living Reviews in Computational Astrophysics 3 (1) (2017) 2.  arXiv:1703.01241, doi:10.1007/s41115-017-0002-8.
  * [32] D. Lee, H. Faller, A. Reyes, The piecewise cubic method (pcm) for computational fluid dynamics, J. Comput. Phys. 341 (2016) 230–257.
  * [33] S. Li, A fourth-order divergence-free method for MHD flows, Journal of Computational Physics 229 (20) (2010) 7893–7910.  doi:10.1016/j.jcp.2010.06.044.
  * [34] J. Núñez-de la Rosa, C.-D. Munz, XTROEM-FV: a new code for computational astrophysics based on very high order finite-volume methods - I. Magnetohydrodynamics, MNRAS455 (4) (2016) 3458–3479.  doi:10.1093/mnras/stv2531.
  * [35] J. Núñez-de la Rosa, C.-D. Munz, xtroem-fv: a new code for computational astrophysics based on very high order finite-volume methods – II. Relativistic hydro- and magnetohydrodynamics, Monthly Notices of the Royal Astronomical Society 460 (1) (2016) 535–559.  arXiv:https://academic.oup.com/mnras/article-pdf/460/1/535/8116559/stw999.pdf, doi:10.1093/mnras/stw999.
URL https://doi.org/10.1093/mnras/stw999

  * [36] K. G. Felker, J. M. Stone, A fourth-order accurate finite volume method for ideal MHD via upwind constrained transport, Journal of Computational Physics 375 (2018) 1365–1400.  arXiv:1711.07439, doi:10.1016/j.jcp.2018.08.025.
  * [37] Y. Matsumoto, Y. Asahina, Y. Kudoh, T. Kawashima, J. Matsumoto, H. R. Takahashi, T. Minoshima, S. Zenitani, T. Miyoshi, R. Matsumoto, Magnetohydrodynamic simulation code CANS+: Assessments and applications, PASJ71 (4) (2019) 83.  arXiv:1611.01775, doi:10.1093/pasj/psz064.
  * [38] P. S. Verma, J.-M. Teissier, O. Henze, W.-C. Müller, Fourth-order accurate finite-volume CWENO scheme for astrophysical MHD problems, MNRAS482 (1) (2019) 416–437.  doi:10.1093/mnras/sty2641.
  * [39] A. Suresh, H. T. Huynh, Accurate Monotonicity-Preserving Schemes with Runge Kutta Time Stepping, Journal of Computational Physics 136 (1997) 83–99.  doi:10.1006/jcph.1997.5745.
  * [40] G.-S. Jiang, C.-W. Shu, Efficient Implementation of Weighted ENO Schemes, Journal of Computational Physics 126 (1) (1996) 202–228.  doi:10.1006/jcph.1996.0130.
  * [41] R. Borges, M. Carmona, B. Costa, W. S. Don, An improved weighted essentially non-oscillatory scheme for hyperbolic conservation laws, Journal of Computational Physics 227 (6) (2008) 3191–3211.  doi:10.1016/j.jcp.2007.11.038.
  * [42] M. Castro, B. Costa, W. S. Don, High order weighted essentially non-oscillatory WENO-Z schemes for hyperbolic conservation laws, Journal of Computational Physics 230 (5) (2011) 1766–1792.  doi:10.1016/j.jcp.2010.11.028.
  * [43] L. Isherwood, Z. J. Grant, S. Gottlieb, Downwinding for Preserving Strong Stability in Explicit Integrating Factor Runge–Kutta Methods, arXiv e-prints (2018) arXiv:1810.04800arXiv:1810.04800.
  * [44] R. J. Spiteri, S. J. Ruuth, A new class of optimal high-order strong-stability-preserving time discretization methods, SIAM Journal on Numerical Analysis 40 (2) (2003) 469–491.
  * [45] C. R. Evans, J. F. Hawley, Simulation of magnetohydrodynamic flows - A constrained transport method, ApJ332 (1988) 659–677.  doi:10.1086/166684.
  * [46] D. Ryu, F. Miniati, T. W. Jones, A. Frank, A Divergence-free Upwind Code for Multidimensional Magnetohydrodynamic Flows, ApJ509 (1) (1998) 244–255.  arXiv:astro-ph/9807228, doi:10.1086/306481.
  * [47] P. Londrillo, L. del Zanna, On the divergence-free condition in Godunov-type schemes for ideal magnetohydrodynamics: the upwind constrained transport method, Journal of Computational Physics 195 (2004) 17–48.  arXiv:astro-ph/0310183, doi:10.1016/j.jcp.2003.09.016.
  * [48] J. U. Brackbill, D. C. Barnes, The Effect of Nonzero $\nabla$ · B on the numerical solution of the magnetohydrodynamic equations, Journal of Computational Physics 35 (3) (1980) 426–430.  doi:10.1016/0021-9991(80)90079-0.
  * [49] K. Yee, Numerical solution of inital boundary value problems involving maxwell’s equations in isotropic media, IEEE Transactions on Antennas and Propagation 14 (3) (1966) 302–307.  doi:10.1109/TAP.1966.1138693.
  * [50] S. H. Brecht, J. Lyon, J. A. Fedder, K. Hain, A simulation study of east-west IMF effects on the magnetosphere, Geophys. Res. Lett.8 (4) (1981) 397–400.  doi:10.1029/GL008i004p00397.
  * [51] V. A. Titarev, E. F. Toro, Finite-volume WENO schemes for three-dimensional conservation laws, Journal of Computational Physics 201 (1) (2004) 238–260.  doi:10.1016/j.jcp.2004.05.015.
  * [52] N. K. Yamaleev, M. H. Carpenter, Third-order Energy Stable WENO scheme, Journal of Computational Physics 228 (2009) 3025–3047.  doi:10.1016/j.jcp.2009.01.011.
  * [53] P. L. Roe, Approximate Riemann Solvers, Parameter Vectors, and Difference Schemes, Journal of Computational Physics 43 (2) (1981) 357–372.  doi:10.1016/0021-9991(81)90128-5.
  * [54] A. Harten, P. Lax, B. Leer, On upstream differencing and godunov-type schemes for hyperbolic conservation laws, SIAM Review 25 (1) (1983) 35–61.  arXiv:https://doi.org/10.1137/1025002, doi:10.1137/1025002.
URL https://doi.org/10.1137/1025002

  * [55] E. F. Toro, M. Spruce, W. Speares, Restoration of the contact surface in the HLL-Riemann solver, Shock Waves 4 (1994) 25–34.  doi:10.1007/BF01414629.
  * [56] S. Li, An HLLC Riemann solver for magneto-hydrodynamics, Journal of Computational Physics 203 (1) (2005) 344–357.  doi:10.1016/j.jcp.2004.08.020.
  * [57] T. Miyoshi, K. Kusano, A multi-state HLL approximate Riemann solver for ideal magnetohydrodynamics, Journal of Computational Physics 208 (1) (2005) 315–344.  doi:10.1016/j.jcp.2005.02.017.
  * [58] A. Mignone, M. Ugliano, G. Bodo, A five-wave Harten-Lax-van Leer Riemann solver for relativistic magnetohydrodynamics, MNRAS393 (2009) 1141–1156.  arXiv:0811.1483, doi:10.1111/j.1365-2966.2008.14221.x.
  * [59] E. Toro, Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction, 2009.  doi:10.1007/b79761.
  * [60] E. F. Toro, V. A. Titarev, MUSTA fluxes for systems of conservation laws, Journal of Computational Physics 216 (2) (2006) 403–429.  doi:10.1016/j.jcp.2005.12.012.
  * [61] C.-W. Shu, S. Osher, Efficient Implementation of Essentially Non-oscillatory Shock-Capturing Schemes, Journal of Computational Physics 77 (2) (1988) 439–471.  doi:10.1016/0021-9991(88)90177-5.
  * [62] S. Gottlieb, C.-W. Shu, E. Tadmor, Strong Stability-Preserving High-Order Time Discretization Methods, SIAM Review 43 (1) (2001) 89–112.  doi:10.1137/S003614450036757X.
  * [63] S. Gottlieb, D. Ketcheson, C. Shu, High order strong stability preserving time discretizations, Journal of Scientific Computing 38 (3) (2009) 251–289.  doi:10.1007/s10915-008-9239-z.
  * [64] S. Gottlieb, D. Ketcheson, C. Shu, Strong Stability Preserving Runge-Kutta and Multistep Time Discretizations, World Scientific, 2011\.
URL https://books.google.it/books?id=MHmAINTBIkQC

  * [65] D. I. Ketcheson, Highly efficient strong stability-preserving runge–kutta methods with low-storage implementations, SIAM Journal on Scientific Computing 30 (4) (2008) 2113–2136.  arXiv:https://doi.org/10.1137/07070485X, doi:10.1137/07070485X.
URL https://doi.org/10.1137/07070485X

  * [66] Numerical Solution of the Euler Equations by Finite Volume Methods Using Runge Kutta Time Stepping Schemes, American Institute of Aeronautics and Astronautics, Palo Alto,CA,U.S.A., 1981.  doi:10.2514/6.1981-1259.
  * [67] M. L. B. de Oliveira, V. A. Pires, Analysis of a high-order finite difference detector for discontinuities, International Journal of Computer Mathematics 94 (4) (2017) 676–689.  arXiv:https://doi.org/10.1080/00207160.2015.1124100, doi:10.1080/00207160.2015.1124100.
URL https://doi.org/10.1080/00207160.2015.1124100

  * [68] A. Mignone, G. Bodo, S. Massaglia, T. Matsakos, O. Tesileanu, C. Zanni, A. Ferrari, PLUTO: A Numerical Code for Computational Astrophysics, ApJS170 (2007) 228–242.  arXiv:astro-ph/0701854, doi:10.1086/513316.
  * [69] A. Mignone, C. Zanni, P. Tzeferacos, B. van Straalen, P. Colella, G. Bodo, The PLUTO Code for Adaptive Mesh Computations in Astrophysical Fluid Dynamics, ApJS198 (2012) 7.  arXiv:1110.0740, doi:10.1088/0067-0049/198/1/7.
  * [70] D. Ryu, T. W. Jones, A. Frank, Numerical Magnetohydrodynamics in Astrophysics: Algorithm and Tests for Multidimensional Flow, ApJ452 (1995) 785\.  arXiv:astro-ph/9505073, doi:10.1086/176347.
  * [71] T. A. Gardiner, J. M. Stone, An unsplit Godunov method for ideal MHD via constrained transport in three dimensions, Journal of Computational Physics 227 (8) (2008) 4123–4141.  arXiv:0712.2634, doi:10.1016/j.jcp.2007.12.017.
  * [72] D. Balsara, Total Variation Diminishing Scheme for Relativistic Magnetohydrodynamics, ApJS132 (1) (2001) 83–101.  doi:10.1086/318941.
  * [73] L. Antón, J. A. Miralles, J. M. Martí, J. M. Ibáñez, M. A. Aloy, P. Mimica, Relativistic Magnetohydrodynamics: Renormalized Eigenvectors and Full Wave Decomposition Riemann Solver, ApJS188 (1) (2010) 1–31.  arXiv:0912.4692, doi:10.1088/0067-0049/188/1/1.
  * [74] M. Dumbser, D. S. Balsara, E. F. Toro, C.-D. Munz, A unified framework for the construction of one-step finite volume and discontinuous Galerkin schemes on unstructured meshes, Journal of Computational Physics 227 (18) (2008) 8209–8253.  doi:10.1016/j.jcp.2008.05.025.
  * [75] A. Mignone, A simple and accurate Riemann solver for isothermal MHD, Journal of Computational Physics 225 (2) (2007) 1427–1441.  arXiv:astro-ph/0701798, doi:10.1016/j.jcp.2007.01.033.
  * [76] K. Beckwith, J. M. Stone, A Second-order Godunov Method for Multi-dimensional Relativistic Magnetohydrodynamics, ApJS193 (1) (2011) 6.  arXiv:1101.3573, doi:10.1088/0067-0049/193/1/6.
  * [77] J. F. Hawley, J. M. Stone, MOCCT: A numerical technique for astrophysical MHD, Computer Physics Communications 89 (1-3) (1995) 127–148.  doi:10.1016/0010-4655(95)00190-Q.
  * [78] S. A. Orszag, C. M. Tang, Small-scale structure of two-dimensional magnetohydrodynamic turbulence, Journal of Fluid Mechanics 90 (1979) 129–143.  doi:10.1017/S002211207900210X.
  * [79] C. F. Gammie, J. C. McKinney, G. Tóth, HARM: A Numerical Scheme for General Relativistic Magnetohydrodynamics, ApJ589 (1) (2003) 444–457.  arXiv:astro-ph/0301509, doi:10.1086/374594.
  * [80] M. Dumbser, O. Zanotti, Very high order PNPM schemes on unstructured meshes for the resistive relativistic MHD equations, Journal of Computational Physics 228 (18) (2009) 6991–7006.  arXiv:0903.4832, doi:10.1016/j.jcp.2009.06.009.
  * [81] S. Landi, L. Del Zanna, E. Papini, F. Pucci, M. Velli, Resistive Magnetohydrodynamics Simulations of the Ideal Tearing Mode, ApJ806 (1) (2015) 131.  arXiv:1504.07036, doi:10.1088/0004-637X/806/1/131.
  * [82] E. Papini, S. Landi, L. Del Zanna, Fast Magnetic Reconnection: Secondary Tearing Instability and Role of the Hall Term, ApJ885 (1) (2019) 56.  arXiv:1906.06779, doi:10.3847/1538-4357/ab4352.
  * [83] G. Mattia, A. Mignone, A comparison of approximate non-linear Riemann solvers for Relativistic MHD, MNRAS510 (1) (2022) 481–499.  arXiv:2111.09369, doi:10.1093/mnras/stab3373.
  * [84] U. Ziegler, A central-constrained transport scheme for ideal magnetohydrodynamics, Journal of Computational Physics 196 (2) (2004) 393–416.  doi:10.1016/j.jcp.2003.11.003.
  * [85] F. Miniati, D. F. Martin, Constrained-transport Magnetohydrodynamics with Adaptive Mesh Refinement in CHARM, ApJS195 (1) (2011) 5.  arXiv:1103.1878, doi:10.1088/0067-0049/195/1/5.