Characteristic boundary conditions for magnetohydrodynamic equations

P. Makaremi-Esfarjani, A. Najafi-Yazdi
Mechanical Engineering Department, McGill University, Montréal, QC, Canada

## Abstract

In the present study, a characteristic-based boundary condition scheme is
developed for the compressible magnetohydrodynamic (MHD) equations in the
general curvilinear coordinate system, which is an extension of the
characteristic boundary scheme for the Navier-Stokes equations. The
eigenstructure and the complete set of characteristic waves are derived for
the ideal MHD equations in general curvilinear coordinates $(\xi,\eta,\zeta)$.
The characteristic boundary conditions are derived and implemented in a high-
order MHD solver where the sixth-order compact scheme is used for the spatial
discretization. The fifth-order Weighted Essentially Non-Oscillatory (WENO)
scheme is also employed for the spatial discretization of problems with
discontinuities. In our MHD solver, the fourth-order Runge-Kutta scheme is
utilized for time integration. The characteristic boundary scheme is first
verified for the non-magnetic (i.e., $\mathbf{B}=\textbf{0}$) Sod shock tube
problem. Then, various in-house test cases are designed to examine the derived
MHD characteristic boundary scheme for three different types of boundaries:
non-reflecting inlet and outlet, solid wall, and single characteristic wave
injection. The numerical examples demonstrate the accuracy and robustness of
the MHD characteristic boundary scheme.
Keywords: Magnetohydrodynamics, MHD, Characteristic boundary, Boundary
conditions, Characteristic waves, Non-reflecting inlet/outlet boundary

##  1 Introduction

Imposing appropriate boundary conditions is of crucial importance to the
accuracy and convergence of numerical results. Characteristic boundary
conditions are among the most accurate boundary conditions used in typical
Computational Fluid Dynamic (CFD) problems involving transient compressible
flows. Such boundary conditions are constructed based on identifying the
outgoing and incoming waves from and to the computational domain. The
amplitude of the incoming waves are set to appropriate values to impose the
desired boundary condition. In this scheme, no extrapolation is used since
different boundary conditions can be applied by using characteristic relations
of the crossing waves at the boundary.

Characteristic boundary conditions for the Euler equations were extensively
studied by Kreiss [20], Engquist and Majda [10], and Hirsch [14]. This scheme
was also studied for the Navier-Stokes equations by Poinsot and Lele [24].
Poinsot and Lele showed that the characteristic method has higher precision
and stability in comparison with other classical approaches such as
extrapolation and simplified Riemann invariants [24, 6]. Moreover, Thompson
[29] studied the time-dependent characteristic boundary conditions for a
general hyperbolic system.

Most numerical tools to solve MHD equations involve periodic computational
domains, which considerably limit the application of such numerical
simulations. One of the most practical boundary conditions is the non-
reflecting inlet/outlet that can be readily implemented using the
characteristic boundary scheme. Non-reflecting boundary conditions for
compressible flows were first studied in the work of Bayliss and Turkel [2] by
using the perturbation approach. In their study, a perturbation was expressed
in terms of the characteristic waves in order to impose boundary conditions
such that it would prevent the formation of incoming waves. The study of
Bayliss and Turkel was followed by the work of Engquist and Majda [10, 11],
which developed the characteristic non-reflecting boundary scheme for linear
acoustic and elastic waves. Later, Hedstrom [13] proposed the non-reflecting
boundary scheme for general nonlinear hyperbolic systems as well. Hedstrom
[13] suggested that a non-reflecting boundary condition involves setting the
amplitude of the incoming waves to zero. Jiang et al. [19] further developed
the non-reflecting boundary conditions in the general curvilinear coordinates
for compressible flows. Here, we expand the idea of the non-reflecting
boundary conditions to the ideal MHD equations.

The application of characteristic boundary conditions for MHD problems has
been somewhat limited. Cimino et al. [6] introduced the characteristic
boundary scheme for the one-dimensional MHD equations in the cartesian
coordinate system. They tried to benchmark two types of boundary conditions,
open end and solid wall, for the Brio–Wu shock tube problem by comparing their
results to that obtained from a zeroth-order extrapolated boundary condition
scheme. This comparison showed some discrepancies that they suggested might be
due to nonphysical boundary conditions, warranting further studies [6].

The goal of the present work is to introduce a set of characteristic boundary
conditions for numerical simulation of the compressible MHD equations,
especially on curvilinear grids. Here, the formal derivation of the
characteristic boundary conditions for the numerical simulation of the MHD
equations is examined in detail. The complete eigenstructure of the MHD
equations is calculated and a set of benchmark problems is designed to verify
the accuracy of the implemented boundary conditions.

##  2 Governing ideal MHD equations

The magnetohydrodynamic set of equations is a fluid model of plasma which
develops a macroscopic description of its dynamics. In this model, plasma is
assumed as a continuous medium where hydrodynamic length scales are much
larger than microscopic ones, namely the particle mean free path and the gyro-
radius. Similarly, hydrodynamic time scales are considered to be much greater
than the mean free time and the gyro-period. Using these assumptions, fluid
dynamic equations can be applied to the motion of plasma. These equations
should be coupled with the Maxwell’s equations to close the system of
equations. The governing equations considered throughout this study are ideal
MHD equations, in which all dissipative processes such as viscosity,
electrical resistivity, and thermal conductivity are neglected.

The ideal MHD equations can be written in the general curvilinear coordinate
system ($\xi$,$\eta$,$\zeta$) as [16]

$\frac{\partial\mathbf{Q}}{\partial
t}+\frac{\partial}{\partial\xi}\left(\frac{\mathbf{F}}{J}\right)+\frac{\partial}{\partial\eta}\left(\frac{\mathbf{G}}{J}\right)+\frac{\partial}{\partial\zeta}\left(\frac{\mathbf{H}}{J}\right)=\frac{\mathbf{S}}{J},$
(1)

where $\mathbf{Q}$ is the vector of the conservative variables,

$\mathbf{Q}=\frac{1}{J}\begin{bmatrix}\rho&\rho u&\rho v&\rho
w&B_{x}&B_{y}&B_{z}&E\end{bmatrix}^{\mathrm{T}}.$ (2)

Variables $J$ and $\rho$ denote the Jacobian and density, respectively.
Velocity field elements in the cartesian coordinate system are shown with $u$,
$v$, and $w$. Furthermore, $B_{x}$, $B_{y}$, and $B_{z}$ represent the
magnetic field components in the physical space, and the total energy value,
$E$, is given as

$E=\frac{p}{\gamma-1}+\rho\frac{\textbf{V}\cdot\textbf{V}}{2}+\frac{\textbf{B}\cdot\textbf{B}}{2},$
(3)

where vectors V and B represent velocity and magnetic field components,
respectively.
The flux terms $\mathbf{F}$, $\mathbf{G}$, and $\mathbf{H}$, are defined as

$\mathbf{F}=\begin{bmatrix}\rho U\\\ \rho
uU+\xi_{x}p+\left(\frac{{B_{x}}^{2}}{2}+\frac{{B_{y}}^{2}}{2}+\frac{{B_{z}}^{2}}{2}\right)\xi_{x}-B_{x}B_{\xi}\\\
\rho
vU+\xi_{y}p+\left(\frac{{B_{x}}^{2}}{2}+\frac{{B_{y}}^{2}}{2}+\frac{{B_{z}}^{2}}{2}\right)\xi_{y}-B_{y}B_{\xi}\\\
\rho
wU+\xi_{z}p+\left(\frac{{B_{x}}^{2}}{2}+\frac{{B_{y}}^{2}}{2}+\frac{{B_{z}}^{2}}{2}\right)\xi_{z}-B_{z}B_{\xi}\\\
B_{x}U-uB_{\xi}\\\ B_{y}U-vB_{\xi}\\\ B_{z}U-wB_{\xi}\\\
\left(E+p+\frac{\mathbf{B}\cdot\mathbf{B}}{2}\right)U-B_{\xi}\left(uB_{x}+vB_{y}+wB_{z}\right)\end{bmatrix},$
(4) $\mathbf{G}=\begin{bmatrix}\rho V\\\ \rho
uV+\eta_{x}p+\left(\frac{{B_{x}}^{2}}{2}+\frac{{B_{y}}^{2}}{2}+\frac{{B_{z}}^{2}}{2}\right)\eta_{x}-B_{x}B_{\eta}\\\
\rho
vV+\eta_{y}p+\left(\frac{{B_{x}}^{2}}{2}+\frac{{B_{y}}^{2}}{2}+\frac{{B_{z}}^{2}}{2}\right)\eta_{y}-B_{y}B_{\eta}\\\
\rho
wV+\eta_{z}p+\left(\frac{{B_{x}}^{2}}{2}+\frac{{B_{y}}^{2}}{2}+\frac{{B_{z}}^{2}}{2}\right)\eta_{z}-B_{z}B_{\eta}\\\
B_{x}V-uB_{\eta}\\\ B_{y}V-vB_{\eta}\\\ B_{z}V-wB_{\eta}\\\
\left(E+p+\frac{\mathbf{B}\cdot\mathbf{B}}{2}\right)V-B_{\eta}\left(uB_{x}+vB_{y}+wB_{z}\right)\end{bmatrix},$
(5)

and

$\mathbf{H}=\begin{bmatrix}\rho W\\\ \rho
uW+\zeta_{x}p+\left(\frac{{B_{x}}^{2}}{2}+\frac{{B_{y}}^{2}}{2}+\frac{{B_{z}}^{2}}{2}\right)\zeta_{x}-B_{x}B_{\zeta}\\\
\rho
vW+\zeta_{y}p+\left(\frac{{B_{x}}^{2}}{2}+\frac{{B_{y}}^{2}}{2}+\frac{{B_{z}}^{2}}{2}\right)\zeta_{y}-B_{y}B_{\zeta}\\\
\rho
wW+\zeta_{z}p+\left(\frac{{B_{x}}^{2}}{2}+\frac{{B_{y}}^{2}}{2}+\frac{{B_{z}}^{2}}{2}\right)\zeta_{z}-B_{z}B_{\zeta}\\\
B_{x}W-uB_{\zeta}\\\ B_{y}W-vB_{\zeta}\\\ B_{z}W-wB_{\zeta}\\\
\left(E+p+\frac{\mathbf{B}\cdot\mathbf{B}}{2}\right)W-B_{\zeta}\left(uB_{x}+vB_{y}+wB_{z}\right)\end{bmatrix},$
(6)

with variables $U$, $V$, and $W$ denoting the velocity components in the
computational space,

$U=u\xi_{x}+v\xi_{y}+w\xi_{z},\ \ V=u\eta_{x}+v\eta_{y}+w\eta_{z},\ \
W=u\zeta_{x}+v\zeta_{y}+w\zeta_{z},$ (7)

where $\xi_{x}$, $\xi_{y}$, and $\xi_{z}$ are curvilinear coordinate metrics.
Additionally, $B_{\xi}$, $B_{\eta}$, and $B_{\zeta}$ variables are the
transformed magnetic field components in the computational space. These
variables are defined in the following equations:

$B_{\xi}=B_{x}\xi_{x}+B_{y}\xi_{y}+B_{z}\xi_{z},\ \
B_{\eta}=B_{x}\eta_{x}+B_{y}\eta_{y}+B_{z}\eta_{z},\ \
B_{\zeta}=B_{x}\zeta_{x}+B_{y}\zeta_{y}+B_{z}\zeta_{z}.$ (8)

Finally, the source term flux vector, $\mathbf{S}$, reads as

$\mathbf{S}=\begin{bmatrix}0\\\ -\left(\frac{\partial
B_{\xi}}{\partial\xi}+\frac{\partial B_{\eta}}{\partial\eta}+\frac{\partial
B_{\zeta}}{\partial\zeta}\right)B_{x}\\\ -\left(\frac{\partial
B_{\xi}}{\partial\xi}+\frac{\partial B_{\eta}}{\partial\eta}+\frac{\partial
B_{\zeta}}{\partial\zeta}\right)B_{y}\\\ -\left(\frac{\partial
B_{\xi}}{\partial\xi}+\frac{\partial B_{\eta}}{\partial\eta}+\frac{\partial
B_{\zeta}}{\partial\zeta}\right)B_{z}\\\ -\left(\frac{\partial
B_{\xi}}{\partial\xi}+\frac{\partial B_{\eta}}{\partial\eta}+\frac{\partial
B_{\zeta}}{\partial\zeta}\right)u\\\ -\left(\frac{\partial
B_{\xi}}{\partial\xi}+\frac{\partial B_{\eta}}{\partial\eta}+\frac{\partial
B_{\zeta}}{\partial\zeta}\right)v\\\ -\left(\frac{\partial
B_{\xi}}{\partial\xi}+\frac{\partial B_{\eta}}{\partial\eta}+\frac{\partial
B_{\zeta}}{\partial\zeta}\right)w\\\ -\left(\frac{\partial
B_{\xi}}{\partial\xi}+\frac{\partial B_{\eta}}{\partial\eta}+\frac{\partial
B_{\zeta}}{\partial\zeta}\right)\left(uB_{x}+vB_{y}+wB_{z}\right)\end{bmatrix}.$
(9)

This source term, known as Powell’s source term, prevents the accumulation of
the finite $\nabla\cdot\mathbf{B}$ generated by the numerical simulation at a
fixed grid point by advecting the error away [25]. This method, also known as
the 8-wave formulation, works well for many problems in which
$\nabla\cdot\mathbf{B}$ remains small [30]. However, for problems containing
strong discontinuities, this method can produce incorrect jump conditions and
inaccurate results. Therefore, the common practice is to couple this method
with a supplementary cleaning technique (e.g., projection scheme) to ensure
that the divergence-free constraint is (approximately) satisfied [30]. Since
the majority of our test cases involve smooth flows and the calculated
$\nabla\cdot\mathbf{B}$ is relatively small, we only apply the 8-wave
formulation scheme for the divergence cleaning.

It is worth mentioning that the derived characteristic boundary scheme in this
study is not only limited to the 8-wave formulation since the source term are
neglected while deriving the MHD eigenstructure. Therefore, the proposed
method can be easily extended to other MHD discretizations using the
projection or constrained transport schemes for satisfying the magnetic field
divergence-free condition. This is due to the fact that the calculated
eigenvalues and eigenvectors for the MHD governing equations remain valid.

##  3 Numerical schemes

###  3.1 Spatial discretization schemes

Two spatial discretization methods are used in this study, namely a sixth-
order central compact scheme [21] and a fifth-order Weighted Essentially Non-
Oscillatory (WENO) scheme [27, 18]. The former is a cost-effective, fast, and
accurate finite difference scheme with spectral-like resolution that is widely
used for problems containing smooth flows. However, this scheme is not
suitable for problems with discontinuities, while the WENO scheme can
accurately capture these discontinuities, maintaining a high order of
accuracy, but at a higher computational cost.

####  3.1.1 Central compact scheme

The sixth-order compact scheme used in this work involves solving a
tridiagonal system of equations for a variable’s derivative. This system is
given by

$\alpha\left(\frac{\partial f}{\partial\xi}\right)_{i-1}+\left(\frac{\partial
f}{\partial\xi}\right)_{i}+\alpha\left(\frac{\partial
f}{\partial\xi}\right)_{i+1}=a\frac{f_{i+1}-f_{i-1}}{2\Delta\xi}+b\frac{f_{i+2}-f_{i-2}}{4\Delta\xi},$
(10)

where $\alpha=1/3$, $a=14/9$, and $b=1/9$.

For the boundary nodes, $i=1$, $N$ and $i=2$, $N-1$, the third-order and
fourth-order one-sided compact schemes were used, respectively. This scheme
has sixth-order accuracy for interior nodes and third- and fourth-order
accuracy for boundary nodes. The central compact scheme can maintain a global
fourth-order accuracy. This tridiagonal system is solved by using the
Tridiagonal Matrix Algorithm (TDMA) [7].

####  3.1.2 WENO scheme

A finite difference version of the WENO scheme [27] is implemented in our MHD
solver code. The main idea of the fifth-order WENO scheme is to use a proper
weighted combination of three sub-stencils for calculating the flux values.
This combination is chosen such that small weight values are assigned to the
stencils with discontinuities to prevent any numerical instability [18]. A
finite volume discretization of the general one-dimensional scalar equation,
$u_{t}=-f(u)_{x}$, can be written as

$\frac{\mathrm{d}u_{i}(t)}{\mathrm{d}t}=-\frac{1}{\Delta
x}(\hat{f}_{i+\frac{1}{2}}-\hat{f}_{i-\frac{1}{2}}),$ (11)

where $u_{i}$ shows the value of $u$ at $x_{i}$, and $\hat{f}$ is a convex
flux which should be constructed in a way to prevent any numerical
instabilities or entropy violating conditions. As a result, the upwinding and
flux splitting approaches should be considered while calculating numerical
fluxes [27]. Here, Lax-Friedrichs splitting [27] is used, giving the positive
and negative flux components as $f^{\pm}(u)=\frac{1}{2}(f(u)\pm\alpha u)$,
where $\alpha$ denotes the maximum value of $|f^{{}^{\prime}}(u)|$. After
finding $f^{\pm}(u)$, we can find the desired positive and negative numerical
fluxes according to the WENO scheme as below [18]:

$\hat{f}^{+}_{i+\frac{1}{2}}=\frac{1}{12}\left(-{f}^{+}_{i-1}+7{f}^{+}_{i}+7{f}^{+}_{i+1}-{f}^{+}_{i+2}\right)-\varphi\left(\Delta{f}^{+}_{i-\frac{3}{2}},\Delta{f}^{+}_{i-\frac{1}{2}},\Delta{f}^{+}_{i+\frac{1}{2}},\Delta{f}^{+}_{i+\frac{3}{2}}\right),$
(12)

and

$\hat{f}^{-}_{i+\frac{1}{2}}=\frac{1}{12}\left(-{f}^{-}_{i-1}+7{f}^{-}_{i}+7{f}^{-}_{i+1}-{f}^{-}_{i+2}\right)+\varphi\left(\Delta{f}^{-}_{i+\frac{5}{2}},\Delta{f}^{-}_{i+\frac{3}{2}},\Delta{f}^{-}_{i+\frac{1}{2}},\Delta{f}^{-}_{i-\frac{1}{2}}\right),$
(13)

where, $\Delta f^{\pm}_{k+\frac{1}{2}}=f^{\pm}_{k+1}-f^{\pm}_{k}$ and
$\varphi$ is a nonlinear function given by [18]

$\varphi(m,n,p,q)=\frac{1}{3}w_{0}\left(m-2n+p\right)+\frac{1}{6}\left(w_{2}-\frac{1}{2}\right)\left(n-2p+q\right).$
(14)

Variables $w_{0}$ and $w_{2}$ are weight values for the first and third
stencils. Jiang and Shu [17] suggested the following weights:

$w_{0}=\frac{\alpha_{0}}{\alpha_{0}+\alpha_{1}+\alpha_{2}},\ \
w_{2}=\frac{\alpha_{2}}{\alpha_{0}+\alpha_{1}+\alpha_{2}},$ (15)

where

$\alpha_{0}=\frac{1}{(\epsilon+\beta_{0})^{2}},\ \
\alpha_{1}=\frac{6}{(\epsilon+\beta_{1})^{2}},\ \
\alpha_{2}=\frac{3}{(\epsilon+\beta_{2})^{2}}.$ (16)

Here parameter $\epsilon$ is used to prevent the denominators to become zero.
The value of $\epsilon$ is typically chosen to be between $10^{-5}$ and
$10^{-7}$. We used a value of $\epsilon=10^{-6}$. The smooth indicator
parameters, denoted by $\beta$, are given as

$\beta_{0}=13(m-n)^{2}+3(m-3n)^{2},\ \ \beta_{1}=13(n-p)^{2}+3(n+p)^{2},\ \
\beta_{2}=13(p-q)^{2}+3(3p-q)^{2}.$ (17)

Combining Eqs. (12) and (13) results in

$\hat{f}_{i+\frac{1}{2}}=\frac{1}{12}\left(-f_{i-1}+7{f}_{i}+7{f}_{i+1}-{f}_{i+2}\right)-\varphi\left(\Delta{f}^{+}_{i-\frac{3}{2}},\Delta{f}^{+}_{i-\frac{1}{2}},\Delta{f}^{+}_{i+\frac{1}{2}},\Delta{f}^{+}_{i+\frac{3}{2}}\right)+\varphi\left(\Delta{f}^{-}_{i+\frac{5}{2}},\Delta{f}^{-}_{i+\frac{3}{2}},\Delta{f}^{-}_{i+\frac{1}{2}},\Delta{f}^{-}_{i-\frac{1}{2}}\right).$
(18)

Similarly, $\hat{f}_{i-\frac{1}{2}}$ can be calculated, and by substituting
$\hat{f}_{i\pm\frac{1}{2}}$ values in Eq. (11), we can calculate the
derivative of $u$ at $x_{i}$.

There are several ways to generalize the scalar WENO scheme to the system of
equations [27]. The most straightforward approach, implemented in this study,
is a component-wise finite difference method [27]. The component by component
version of the WENO scheme is simple, cost-effective, and reasonably robust
for many problems. In this method, the flux splitting procedure is applied
separately to each component of $u$ to construct the numerical flux values
${\hat{f}}_{i+\frac{1}{2}}$. Afterward, the spatial derivative is calculated
using Eq. (11). Furthermore, the same procedure should be applied in each
direction for multi-dimensional problems.

###  3.2 Temporal integration scheme

Consider a system of differential equations given by
$\frac{\partial\mathbf{U}}{\partial t}=\mathbf{F}(\mathbf{U},t)$, where
$\mathbf{U}$ is the vector of conservative variables and $\mathbf{F}$ is a
right-hand-side operator. For a known solution, $\mathbf{U}^{n}$, the low-
storage format of the fourth-order Runge-Kutta scheme [22] is implemented in
our in-house code to approximate the solution at the next time step,
$\mathbf{U}^{n+1}$. This method uses the following four steps to approximate
the solution at the time step $n+1$:

  1. 1.

The first step (Euler Predictor):

$\displaystyle\mathbf{U}^{{}^{\prime}}=\mathbf{F}(\mathbf{U}^{n},t)$ (19a)
$\displaystyle\mathbf{Q}_{1}=\mathbf{U}^{n}+\frac{\Delta
t}{2}\mathbf{U}^{{}^{\prime}}$ (19b)
$\displaystyle\mathbf{Q}_{2}=\mathbf{U}^{n}+\frac{\Delta
t}{6}\mathbf{U}^{{}^{\prime}}$ (19c)

  2. 2.

The second step (Euler Corrector):

$\displaystyle\mathbf{U}^{{}^{\prime}}=\mathbf{F}\left(\mathbf{Q}_{1},t+\frac{\Delta
t}{2}\right)$ (20a) $\displaystyle\mathbf{Q}_{1}=\mathbf{U}^{n}+\frac{\Delta
t}{2}\mathbf{U}^{{}^{\prime}}$ (20b)
$\displaystyle\mathbf{Q}_{2}=\mathbf{Q}_{2}+\frac{\Delta
t}{3}\mathbf{U}^{{}^{\prime}}$ (20c)

  3. 3.

The third step (Leapfrog Predictor):

$\displaystyle\mathbf{U}^{{}^{\prime}}=\mathbf{F}\left(\mathbf{Q}_{1},t+\frac{\Delta
t}{2}\right)$ (21a) $\displaystyle\mathbf{Q}_{1}=\mathbf{U}^{n}+\Delta
t\mathbf{U}^{{}^{\prime}}$ (21b)
$\displaystyle\mathbf{Q}_{2}=\mathbf{Q}_{2}+\frac{\Delta
t}{3}\mathbf{U}^{{}^{\prime}}$ (21c)

  4. 4.

The fourth step (Milne Corrector) :

$\displaystyle\mathbf{U}^{{}^{\prime}}=\mathbf{F}\left(\mathbf{Q}_{1},t+\Delta
t\right)$ (22a) $\displaystyle\mathbf{U}^{n+1}=\mathbf{Q}_{2}+\frac{\Delta
t}{6}\mathbf{U}^{{}^{\prime}}$ (22b)

The convergence study of the implemented MHD solver is provided in Appendix A
for interested readers.

##  4 Characteristic waves and boundary conditions

In this section, the procedure of implementing the characteristic boundary
scheme for a one-dimensional system of governing equations is presented. The
same methodology will be extended to the MHD equations in three-dimensional
curvilinear coordinates using the complete derived set of MHD eigenstructure
provided in Sec. 4.2. In Sec. 4.3, we delineate how different boundary
conditions, namely non-reflecting inlet/outlet, solid wall, and constant
pressure, can be imposed by employing a characteristic boundary scheme.

###  4.1 Characteristic boundary scheme

The characteristic boundary scheme specifies the desired boundary conditions
by identifying and manipulating the different waves crossing the computational
boundaries [23]. The information from the outgoing waves are used to satisfy
the desired boundary conditions by imposing proper constraints on the incoming
waves [6].

Consider the conservative form of the one-dimensional system of governing
equations in the cartesian coordinate system

$\frac{\partial\mathbf{U}}{\partial
t}+\frac{\partial\mathbf{F}(\mathbf{U})}{\partial x}=\textbf{0},$ (23)

where $\mathbf{U}$ is the vector of conservative variables and $\mathbf{F}$ is
the flux vector. Eq. (23) describes the relationship between the temporal rate
of the conservative field over an infinitesimal control volume and the flux of
that field crossing the control volume boundaries [29]. Eq. (23) can also be
written in primitive form as [29]:

$\frac{\partial\mathbf{V}}{\partial
t}+\mathbf{A}\frac{\partial\mathbf{V}}{\partial x}=\textbf{0},$ (24)

where $\mathbf{V}$ is the vector of primitive variables and $\mathbf{A}$ is
known as the Jacobian matrix of the primitive variables in the $x$-direction.
By defining $\mathbf{P}=\partial\mathbf{U}/\partial\mathbf{V}$ and
$\mathbf{Q}=\partial\mathbf{F}/\partial\mathbf{V}$, we can rewrite Eq. (23) as
[24, 29, 13]:

$\frac{\partial\mathbf{U}}{\partial
t}+\mathbf{Q}\frac{\partial\mathbf{V}}{\partial{x}}=\textbf{0}.$ (25)

Multiplying Eq. (25) by $\mathbf{P}^{-1}$ will result in Eq. (24) in which
$\mathbf{P}^{-1}\mathbf{Q}=\mathbf{A}$. Therefore, for converting the
conservative form of the governing equations to the primitive form, we should
multiply Eq. (23) by $\mathbf{P}^{-1}$. Matrix $\mathbf{A}$ is a
diagonalizable matrix and can be written in the form
$\mathbf{A}=\mathbf{R}\mathbf{\Lambda}\mathbf{L}$, where $\mathbf{R}$ and
$\mathbf{L}$ are the right and left eigenvectors, and $\mathbf{\Lambda}$ is
the diagonal matrix of the corresponding eigenvalues. Using this definition,
we can express Eq. (23) as

$\frac{\partial\mathbf{V}}{\partial
t}+\mathbf{R}\mathbf{\Lambda}\mathbf{L}\frac{\partial\mathbf{V}}{\partial
x}=\textbf{0}.$ (26)

By multiplying both sides of Eq. (26) by the left eigenvector and noting that
$\mathbf{LR}=\mathbf{I}$ (identity matrix), we obtain

$\mathbf{L}\frac{\partial\mathbf{V}}{\partial
t}+\Lambda\mathbf{L}\frac{\partial\mathbf{V}}{\partial x}=\textbf{0},\ \
\text{or}\ \ \mathbf{l}_{i}\frac{\partial\mathbf{V}}{\partial
t}+\lambda_{i}\mathbf{l}_{i}\frac{\partial\mathbf{V}}{\partial x}=\textbf{0},$
(27)

where $\mathbf{l}_{i}$ denotes the $i$-th row of the left eigenvector and
$\lambda_{i}$ is its corresponding eigenvalue. Equation (27) shows that by
projecting the primitive variables in the eigenvector space, we can decouple
the governing equations which can be treated individually. For simplicity, we
shall denote $\lambda_{i}\mathbf{l}_{i}\frac{\partial\mathbf{V}}{\partial x}$
by operators $\mathcal{L}_{i}$ hereafter. In order to express the derivative
of the primitive variables as a function of $\mathcal{L}_{i}$ operators, we
multiply Eq. (27) by the right eigenvectors. This operation will result in

$\frac{\partial\mathbf{V}}{\partial t}+\mathbf{R}\mathcal{L}=\textbf{0},\ \
\text{or}\ \ \frac{\partial\mathbf{V}}{\partial
t}+\mathbf{r}_{i}\mathcal{L}_{i}=\textbf{0},$ (28)

where $\mathbf{r}_{i}$ denotes the columns of the right eigenvector. At this
step, the time evolution of the primitive variables at the boundaries is
derived as a function of $\mathcal{L}_{i}$ operators. Using this set of
equations, we can impose the desired boundary conditions by modifying
$\mathcal{L}_{i}$ operators accordingly, i.e.,
$\mathbf{L}\,\partial\mathbf{V}/\partial
t+\mathcal{L}_{\text{modified}}=\textbf{0}$. We can convert this modified
equation to the conservative form by multiplying it by the right eigenvector
and $\mathbf{P}$, respectively. This operation will result in
$\partial\mathbf{U}/\partial t+\left(\partial\mathbf{F}(\mathbf{U})/\partial
x\right)_{\text{modified}}=\textbf{0}$, which implements the desired boundary
condition by modifying the right-hand-side operator of the Runge-Kutta scheme
(see section 3.2). Figure 1 shows the procedure for implementing the
characteristic boundary scheme for the system of governing equations.

Fig. 1: A flowchart for implementing the characteristic boundary scheme for
the general system of governing equations.

It is worth noting that the described procedure of implementing the
characteristic boundary condition is stable in discretization. This is due to
the fact that the characteristic boundary treatment is done after the flux
discretization. Therefore, it would not create any instability in the
numerical discretization.

###  4.2 MHD eigenstructure in the general curvilinear coordinate system

Consider a three-dimensional set of MHD equations written in the general form
of

$\frac{\partial\mathbf{V}}{\partial
t}+\mathbf{A}_{\xi}\frac{\partial\mathbf{V}}{\partial\xi}+\mathbf{A}_{\eta}\frac{\partial\mathbf{V}}{\partial\eta}+\mathbf{A}_{\zeta}\frac{\partial\mathbf{V}}{\partial\zeta}=\textbf{0}.$
(29)

Matrices $\mathbf{A}_{\xi}$, $\mathbf{A}_{\eta}$, and $\mathbf{A}_{\zeta}$ are
Jacobians of the primitive variables in $\xi$, $\eta$, and $\zeta$ directions
in the general curvilinear coordinate system.

The matrix $\mathbf{A}_{\xi}$ is defined as

$\mathbf{A}_{\xi}=\begin{bmatrix}U&\rho\xi_{x}&\rho\xi_{y}&\rho\xi_{z}&0&0&0&0\\\
0&U&0&0&-\frac{B_{y}}{\rho}\xi_{y}-\frac{B_{z}}{\rho}\xi_{z}&\frac{B_{y}}{\rho}\xi_{x}&\frac{B_{z}}{\rho}\xi_{x}&\frac{\xi_{x}}{\rho}\\\
0&0&U&0&\frac{B_{x}}{\rho}\xi_{y}&-\frac{B_{x}}{\rho}\xi_{x}-\frac{B_{z}}{\rho}\xi_{z}&\frac{B_{z}}{\rho}\xi_{y}&\frac{\xi_{y}}{\rho}\\\
0&0&0&U&\frac{B_{x}}{\rho}\xi_{z}&\frac{B_{y}}{\rho}\xi_{z}&-\frac{B_{x}}{\rho}\xi_{x}-\frac{B_{y}}{\rho}\xi_{y}&\frac{\xi_{z}}{\rho}\\\
0&-B_{y}\xi_{y}-B_{z}\xi_{z}&B_{x}\xi_{y}&B_{x}\xi_{z}&U&0&0&0\\\
0&B_{y}\xi_{x}&-B_{x}\xi_{x}-B_{z}\xi_{z}&B_{y}\xi_{z}&0&U&0&0\\\
0&B_{z}\xi_{x}&B_{z}\xi_{y}&-B_{x}\xi_{x}-B_{y}\xi_{y}&0&0&U&0\\\ 0&\gamma
p\xi_{x}&\gamma p\xi_{y}&\gamma p\xi_{z}&0&0&0&U\end{bmatrix},$ (30)

where $U=u\xi_{x}+v\xi_{y}+w\xi_{z}$. Matrices $\mathbf{A}_{\eta}$ and
$\mathbf{A}_{\zeta}$ can be defined akin to Eq. (30) by replacing the variable
$\xi$ with $\eta$ and $\zeta$, respectively. The eigenvalues of the matrix
$\mathbf{A}_{\xi}$ are

where $c_{\mathrm{f}}$ and $c_{\mathrm{s}}$ are the fast and slow magnetosonic
waves derived in the general curvilinear coordinate system. These variables
are defined as

${c_{\mathrm{f}}}^{2}=\frac{1}{2}\left[\left({{a}}^{2}+\left|{{b}}^{2}\right|\right)+\sqrt{\left[{\left({a}^{2}+\left|{{b}}^{2}\right|\right)^{2}-4{{a}}^{2}{{b_{\xi}}}^{2}}\right]}\right],\
\
{c_{\mathrm{s}}}^{2}=\frac{1}{2}\left[\left({{a}}^{2}+\left|{{b}}^{2}\right|\right)-\sqrt{\left[{\left({a}^{2}+\left|{{b}}^{2}\right|\right)^{2}-4{{a}}^{2}{{b_{\xi}}}^{2}}\right]}\right],$
(31)

where

${{a}}^{2}=\frac{\gamma
p}{\rho}\left({\xi_{x}}^{2}+{\xi_{y}}^{2}+{\xi_{z}}^{2}\right),\ \
\left|{{b}}^{2}\right|=\left({b_{x}}^{2}+{b_{y}}^{2}+{b_{z}}^{2}\right)\left({\xi_{x}}^{2}+{\xi_{y}}^{2}+{\xi_{z}}^{2}\right),\
\ {b_{\xi}}=b_{x}\xi_{x}+b_{y}\xi_{y}+b_{z}\xi_{z},\ \
\mathbf{b}=\mathbf{B}/\sqrt{\rho}.$ (32)

Variable ${c_{\mathrm{a}}}$ is the Alfvén wave speed in the curvilinear
coordinates, which is given by

${c_{\mathrm{a}}}={{b_{\xi}}}.$ (33)

Eigenvalues ${\lambda_{1}}$ and ${\lambda_{8}}$ are the propagation speed of
the fast magnetosonic characteristic waves, while eigenvalues ${\lambda_{3}}$
and ${\lambda_{6}}$ are related to the slow magnetosonic waves. Eigenvalues
${\lambda_{2}}$ and ${\lambda_{7}}$ correspond to the Alfvén characteristic
waves. There are two equal eigenvalues ${\lambda_{4}}={\lambda_{5}}=u$ in Eq.
(4.2). The former is the propagation speed of the entropy wave, whereas the
latter corresponds to a fictitious wave, convecting $\nabla\cdot\mathbf{B}$
error.

Magnetosonic waves (slow and fast) are analogous to the acoustic ones that are
modified by the presence of the magnetic field. These waves are longitudinal
and compressible. The restoring force consists of the magnetic pressure,
$B_{i}B_{i}$, and the plasma pressure, $p$. In the case of the slow
magnetosonic mode, density and the absolute value of the magnetic field
perturbations are out of phase, while they are in phase for the fast
magnetosonic modes. Alfvén waves propagate in a direction parallel to the
magnetic field. The magnetic tension force, $B_{i}B_{j}$, is known as a
restoring force in this case. Alfvén waves do not affect the density field and
hence are considered incompressible. The right eigenvector of
$\mathbf{A}_{\xi}$ is

$\displaystyle{r_{1}}=\begin{array}[]{ccccc}\bigl{[}\rho{\alpha_{\mathrm{f}}}&\beta_{x}c_{\mathrm{s}}\alpha_{\mathrm{s}}S-\frac{\xi_{x}c_{\mathrm{f}}\alpha_{\mathrm{s}}a}{b_{\perp}{h_{\xi}}^{2}}&\beta_{y}c_{\mathrm{s}}\alpha_{\mathrm{s}}S-\frac{\xi_{y}c_{\mathrm{f}}\alpha_{\mathrm{s}}a}{b_{\perp}{h_{\xi}}^{2}}&\beta_{z}c_{\mathrm{s}}\alpha_{\mathrm{s}}S-\frac{\xi_{z}c_{\mathrm{f}}\alpha_{\mathrm{s}}a}{b_{\perp}{h_{\xi}}^{2}}&\beta_{x}a\alpha_{\mathrm{s}}\sqrt{\rho}-\frac{\xi_{x}\beta_{\xi}a\alpha_{\mathrm{s}}\sqrt{\rho}}{{h_{\xi}}^{2}}\\\
&\beta_{y}a\alpha_{\mathrm{s}}\sqrt{\rho}-\frac{\xi_{y}\beta_{\xi}a\alpha_{\mathrm{s}}\sqrt{\rho}}{{h_{\xi}}^{2}}&\beta_{z}a\alpha_{\mathrm{s}}\sqrt{\rho}-\frac{\xi_{z}\beta_{\xi}a\alpha_{\mathrm{s}}\sqrt{\rho}}{{h_{\xi}}^{2}}&\frac{\rho
a^{2}\alpha_{\mathrm{f}}}{{h_{\xi}}^{2}}\bigr{]},\end{array}$ (34)
$\displaystyle{r_{2}}=\begin{array}[]{ccccc}\bigl{[}0&\beta_{z}\xi_{y}S-\beta_{y}\xi_{z}S&\beta_{x}\xi_{z}S-\beta_{z}\xi_{x}S&\beta_{y}\xi_{x}S-\beta_{x}\xi_{y}S&\sqrt{\rho}\beta_{z}\xi_{y}-\sqrt{\rho}\beta_{y}\xi_{z}\\\
&\sqrt{\rho}\beta_{x}\xi_{z}-\sqrt{\rho}\beta_{z}\xi_{x}&\sqrt{\rho}\beta_{y}\xi_{x}-\sqrt{\rho}\beta_{x}\xi_{y}&0\bigr{]},\end{array}$
$\displaystyle{r_{3}}=\begin{array}[]{ccccc}\bigl{[}\rho{\alpha_{\mathrm{s}}}&-\beta_{x}c_{\mathrm{f}}\alpha_{\mathrm{f}}S+\frac{\xi_{x}c_{\mathrm{s}}\alpha_{\mathrm{f}}a}{b_{\perp}{h_{\xi}}^{2}}&-\beta_{y}c_{\mathrm{f}}\alpha_{\mathrm{f}}S+\frac{\xi_{y}c_{\mathrm{s}}\alpha_{\mathrm{f}}a}{b_{\perp}{h_{\xi}}^{2}}&-\beta_{z}c_{\mathrm{f}}\alpha_{\mathrm{f}}S+\frac{\xi_{z}c_{\mathrm{s}}\alpha_{\mathrm{f}}a}{b_{\perp}{h_{\xi}}^{2}}&-\beta_{x}a\alpha_{\mathrm{f}}\sqrt{\rho}+\frac{\xi_{x}\beta_{\xi}a\alpha_{\mathrm{f}}\sqrt{\rho}}{{h_{\xi}}^{2}}\\\
&-\beta_{y}a\alpha_{\mathrm{f}}\sqrt{\rho}+\frac{\xi_{y}\beta_{\xi}a\alpha_{\mathrm{f}}\sqrt{\rho}}{{h_{\xi}}^{2}}&-\beta_{z}a\alpha_{\mathrm{f}}\sqrt{\rho}+\frac{\xi_{z}\beta_{\xi}a\alpha_{\mathrm{f}}\sqrt{\rho}}{{h_{\xi}}^{2}}&\frac{\rho
a^{2}\alpha_{\mathrm{s}}}{{h_{\xi}}^{2}}\bigr{]},\end{array}$
$\displaystyle{r_{4}}=\begin{bmatrix}1&0&0&0&0&0&0&0\end{bmatrix},$
$\displaystyle{r_{5}}=\begin{bmatrix}0&0&0&0&\xi_{x}&\xi_{y}&\xi_{z}&0\end{bmatrix},$
$\displaystyle{r_{6}}=\begin{array}[]{ccccc}\bigl{[}\rho{\alpha_{\mathrm{s}}}&\beta_{x}c_{\mathrm{f}}\alpha_{\mathrm{f}}S-\frac{\xi_{x}c_{\mathrm{s}}\alpha_{\mathrm{f}}a}{b_{\perp}{h_{\xi}}^{2}}&\beta_{y}c_{\mathrm{f}}\alpha_{\mathrm{f}}S-\frac{\xi_{y}c_{\mathrm{s}}\alpha_{\mathrm{f}}a}{b_{\perp}{h_{\xi}}^{2}}&\beta_{z}c_{\mathrm{f}}\alpha_{\mathrm{f}}S-\frac{\xi_{z}c_{\mathrm{s}}\alpha_{\mathrm{f}}a}{b_{\perp}{h_{\xi}}^{2}}&-\beta_{x}a\alpha_{\mathrm{f}}\sqrt{\rho}+\frac{\xi_{x}\beta_{\xi}a\alpha_{\mathrm{f}}\sqrt{\rho}}{{h_{\xi}}^{2}}\\\
&-\beta_{y}a\alpha_{\mathrm{f}}\sqrt{\rho}+\frac{\xi_{y}\beta_{\xi}a\alpha_{\mathrm{f}}\sqrt{\rho}}{{h_{\xi}}^{2}}&-\beta_{z}a\alpha_{\mathrm{f}}\sqrt{\rho}+\frac{\xi_{z}\beta_{\xi}a\alpha_{\mathrm{f}}\sqrt{\rho}}{{h_{\xi}}^{2}}&\frac{\rho
a^{2}\alpha_{\mathrm{s}}}{{h_{\xi}}^{2}}\bigr{]},\end{array}$
$\displaystyle{r_{7}}=\begin{array}[]{ccccc}\bigl{[}0&-\beta_{z}\xi_{y}S+\beta_{y}\xi_{z}S&-\beta_{x}\xi_{z}S+\beta_{z}\xi_{x}S&-\beta_{y}\xi_{x}S+\beta_{x}\xi_{y}S&\sqrt{\rho}\beta_{z}\xi_{y}-\sqrt{\rho}\beta_{y}\xi_{z}\\\
&\sqrt{\rho}\beta_{x}\xi_{z}-\sqrt{\rho}\beta_{z}\xi_{x}&\sqrt{\rho}\beta_{y}\xi_{x}-\sqrt{\rho}\beta_{x}\xi_{y}&0\bigr{]},\end{array}$
$\displaystyle{r_{8}}=\begin{array}[]{ccccc}\bigl{[}\rho{\alpha_{\mathrm{f}}}&-\beta_{x}c_{\mathrm{s}}\alpha_{\mathrm{s}}S+\frac{\xi_{x}c_{\mathrm{f}}\alpha_{\mathrm{s}}a}{b_{\perp}{h_{\xi}}^{2}}&-\beta_{y}c_{\mathrm{s}}\alpha_{\mathrm{s}}S+\frac{\xi_{y}c_{\mathrm{f}}\alpha_{\mathrm{s}}a}{b_{\perp}{h_{\xi}}^{2}}&-\beta_{z}c_{\mathrm{s}}\alpha_{\mathrm{s}}S+\frac{\xi_{z}c_{\mathrm{f}}\alpha_{\mathrm{s}}a}{b_{\perp}{h_{\xi}}^{2}}&\beta_{x}a\alpha_{\mathrm{s}}\sqrt{\rho}-\frac{\xi_{x}\beta_{\xi}a\alpha_{\mathrm{s}}\sqrt{\rho}}{{h_{\xi}}^{2}}\\\
&\beta_{y}a\alpha_{\mathrm{s}}\sqrt{\rho}-\frac{\xi_{y}\beta_{\xi}a\alpha_{\mathrm{s}}\sqrt{\rho}}{{h_{\xi}}^{2}}&\beta_{z}a\alpha_{\mathrm{s}}\sqrt{\rho}-\frac{\xi_{z}\beta_{\xi}a\alpha_{\mathrm{s}}\sqrt{\rho}}{{h_{\xi}}^{2}}&\frac{\rho
a^{2}\alpha_{\mathrm{f}}}{{h_{\xi}}^{2}}\bigr{]},\end{array}$

where ${r_{1}}$ to ${r_{8}}$ are the columns of the right eigenvector
corresponding to the eigenvalues ${\lambda_{1}}$ to ${\lambda_{8}}$,
respectively. Dimensionless variables $\alpha_{\mathrm{s}}$,
$\alpha_{\mathrm{f}}$, $\beta$, and $S$ are defined as:

${{\alpha_{\mathrm{s}}}}^{2}=\frac{{{c_{\mathrm{f}}}}^{2}-{{a}}^{2}}{{{c_{\mathrm{f}}}}^{2}-{{c_{\mathrm{s}}}}^{2}},\
\
{{\alpha_{\mathrm{f}}}}^{2}=\frac{{{a}}^{2}-{{c_{\mathrm{s}}}}^{2}}{{c_{\mathrm{f}}}^{2}-{{c_{\mathrm{s}}}}^{2}},\
\ {\beta_{x(y,z)}}=\frac{{b_{x(y,z)}}}{{b_{\perp}}},\ \
{\beta_{\xi}}=\frac{{b_{\xi}}}{{b_{\perp}}},\ \ {S}=\mathrm{sign}({b_{\xi}}).$
(35)

Moreover, variables $b_{\perp}$ and ${h_{\xi}}^{2}$ are given as

${b_{\perp}}=\sqrt{{{b}}^{2}-{{b_{\xi}}}^{2}},\ \
{h_{\xi}}^{2}={\xi_{x}}^{2}+{\xi_{y}}^{2}+{\xi_{z}}^{2}.$ (36)

Here, parameter ${\alpha_{\mathrm{f}}}$ (${\alpha_{\mathrm{s}}}$) shows how
close fast (slow) magnetosonic waves behave in comparison with the acoustic
waves [26]. The rows of the left eigenvector corresponding to the eigenvalues
${\lambda_{1}}$ to ${\lambda_{8}}$ are calculated as

$\displaystyle{l_{1}}=\begin{array}[]{ccccc}\bigl{[}0&-\frac{c_{\mathrm{f}}\alpha_{\mathrm{s}}\xi_{x}}{2ab_{\perp}}+\frac{\beta_{x}c_{\mathrm{s}}\alpha_{\mathrm{s}}{h_{\xi}}^{2}S}{2a^{2}}&-\frac{c_{\mathrm{f}}\alpha_{\mathrm{s}}\xi_{y}}{2ab_{\perp}}+\frac{\beta_{y}c_{\mathrm{s}}\alpha_{\mathrm{s}}{h_{\xi}}^{2}S}{2a^{2}}&-\frac{c_{\mathrm{f}}\alpha_{\mathrm{s}}\xi_{z}}{2ab_{\perp}}+\frac{\beta_{z}c_{\mathrm{s}}\alpha_{\mathrm{s}}{h_{\xi}}^{2}S}{2a^{2}}&\frac{\beta_{x}\alpha_{\mathrm{s}}{h_{\xi}}^{2}-\beta_{\xi}\alpha_{\mathrm{s}}\xi_{x}}{2a\sqrt{\rho}}\\\
&\frac{\beta_{y}\alpha_{\mathrm{s}}{h_{\xi}}^{2}-\beta_{\xi}\alpha_{\mathrm{s}}\xi_{y}}{2a\sqrt{\rho}}&\frac{\beta_{z}\alpha_{\mathrm{s}}{h_{\xi}}^{2}-\beta_{\xi}\alpha_{\mathrm{s}}\xi_{z}}{2a\sqrt{\rho}}&\frac{\alpha_{\mathrm{f}}{h_{\xi}}^{2}}{2\rho
a^{2}}\bigr{]},\end{array}$ (37)
$\displaystyle{l_{2}}=\begin{array}[]{ccccccc}\bigl{[}0&\frac{\beta_{z}\xi_{y}-\beta_{y}\xi_{z}}{2}&\frac{\beta_{x}\xi_{z}-\beta_{z}\xi_{x}}{2}&\frac{\beta_{y}\xi_{x}-\beta_{x}\xi_{y}}{2}&\frac{\beta_{z}\xi_{y}S-\beta_{y}\xi_{z}S}{2\sqrt{\rho}}&\frac{\beta_{x}\xi_{z}S-\beta_{z}\xi_{x}S}{2\sqrt{\rho}}\\\
&\frac{\beta_{y}\xi_{x}S-\beta_{x}\xi_{y}S}{2\sqrt{\rho}}&0\bigr{]},\end{array}$
$\displaystyle{l_{3}}=\begin{array}[]{ccccc}\bigl{[}0&\frac{c_{\mathrm{s}}\alpha_{\mathrm{f}}\xi_{x}}{2ab_{\perp}}-\frac{\beta_{x}c_{\mathrm{f}}\alpha_{\mathrm{f}}{h_{\xi}}^{2}S}{2a^{2}}&\frac{c_{\mathrm{s}}\alpha_{\mathrm{f}}\xi_{y}}{2ab_{\perp}}-\frac{\beta_{y}c_{\mathrm{f}}\alpha_{\mathrm{f}}{h_{\xi}}^{2}S}{2a^{2}}&\frac{c_{\mathrm{s}}\alpha_{\mathrm{f}}\xi_{z}}{2ab_{\perp}}-\frac{\beta_{z}c_{\mathrm{f}}\alpha_{\mathrm{f}}{h_{\xi}}^{2}S}{2a^{2}}&\frac{-\beta_{x}\alpha_{\mathrm{f}}{h_{\xi}}^{2}+\beta_{\xi}\alpha_{\mathrm{f}}\xi_{x}}{2a\sqrt{\rho}}\\\
&\frac{-\beta_{y}\alpha_{\mathrm{f}}{h_{\xi}}^{2}+\beta_{\xi}\alpha_{\mathrm{f}}\xi_{y}}{2a\sqrt{\rho}}&\frac{-\beta_{z}\alpha_{\mathrm{f}}{h_{\xi}}^{2}+\beta_{\xi}\alpha_{\mathrm{f}}\xi_{z}}{2a\sqrt{\rho}}&\frac{\alpha_{\mathrm{s}}{h_{\xi}}^{2}}{2\rho
a^{2}}\bigr{]},\end{array}$
$\displaystyle{l_{4}}=\begin{bmatrix}1&0&0&0&0&0&0&-\frac{{h_{\xi}}^{2}}{{{a}}^{2}}\end{bmatrix},$
$\displaystyle{l_{5}}=\begin{bmatrix}0&0&0&0&\frac{\xi_{x}}{{h_{\xi}}^{2}}&\frac{\xi_{y}}{{h_{\xi}}^{2}}&\frac{\xi_{z}}{{h_{\xi}}^{2}}&0\end{bmatrix},$
$\displaystyle{l_{6}}=\begin{array}[]{ccccc}\bigl{[}0&-\frac{c_{\mathrm{s}}\alpha_{\mathrm{f}}\xi_{x}}{2ab_{\perp}}+\frac{\beta_{x}c_{\mathrm{f}}\alpha_{\mathrm{f}}{h_{\xi}}^{2}S}{2a^{2}}&-\frac{c_{\mathrm{s}}\alpha_{\mathrm{f}}\xi_{y}}{2ab_{\perp}}+\frac{\beta_{y}c_{\mathrm{f}}\alpha_{\mathrm{f}}{h_{\xi}}^{2}S}{2a^{2}}&-\frac{c_{\mathrm{s}}\alpha_{\mathrm{f}}\xi_{z}}{2ab_{\perp}}+\frac{\beta_{z}c_{\mathrm{f}}\alpha_{\mathrm{f}}{h_{\xi}}^{2}S}{2a^{2}}&\frac{-\beta_{x}\alpha_{\mathrm{f}}{h_{\xi}}^{2}+\beta_{\xi}\alpha_{\mathrm{f}}\xi_{x}}{2a\sqrt{\rho}}\\\
&\frac{-\beta_{y}\alpha_{\mathrm{f}}{h_{\xi}}^{2}+\beta_{\xi}\alpha_{\mathrm{f}}\xi_{y}}{2a\sqrt{\rho}}&\frac{-\beta_{z}\alpha_{\mathrm{f}}{h_{\xi}}^{2}+\beta_{\xi}\alpha_{\mathrm{f}}\xi_{z}}{2a\sqrt{\rho}}&\frac{\alpha_{\mathrm{s}}{h_{\xi}}^{2}}{2\rho
a^{2}}\bigr{]},\end{array}$
$\displaystyle{l_{7}}=\begin{array}[]{ccccccc}\bigl{[}0&\frac{-\beta_{z}\xi_{y}+\beta_{y}\xi_{z}}{2}&\frac{-\beta_{x}\xi_{z}+\beta_{z}\xi_{x}}{2}&\frac{-\beta_{y}\xi_{x}+\beta_{x}\xi_{y}}{2}&\frac{\beta_{z}\xi_{y}S-\beta_{y}\xi_{z}S}{2\sqrt{\rho}}&\frac{\beta_{x}\xi_{z}S-\beta_{z}\xi_{x}S}{2\sqrt{\rho}}\\\
&\frac{\beta_{y}\xi_{x}S-\beta_{x}\xi_{y}S}{2\sqrt{\rho}}&0\bigr{]},\end{array}$
$\displaystyle{l_{8}}=\begin{array}[]{ccccc}\bigl{[}0&\frac{c_{\mathrm{f}}\alpha_{\mathrm{s}}\xi_{x}}{2ab_{\perp}}-\frac{\beta_{x}c_{\mathrm{s}}\alpha_{\mathrm{s}}{h_{\xi}}^{2}S}{2a^{2}}&\frac{c_{\mathrm{f}}\alpha_{\mathrm{s}}\xi_{y}}{2ab_{\perp}}-\frac{\beta_{y}c_{\mathrm{s}}\alpha_{\mathrm{s}}{h_{\xi}}^{2}S}{2a^{2}}&\frac{c_{\mathrm{f}}\alpha_{\mathrm{s}}\xi_{z}}{2ab_{\perp}}-\frac{\beta_{z}c_{\mathrm{s}}\alpha_{\mathrm{s}}{h_{\xi}}^{2}S}{2a^{2}}&\frac{\beta_{x}\alpha_{\mathrm{s}}{h_{\xi}}^{2}-\beta_{\xi}\alpha_{\mathrm{s}}\xi_{x}}{2a\sqrt{\rho}}\\\
&\frac{\beta_{y}\alpha_{\mathrm{s}}{h_{\xi}}^{2}-\beta_{\xi}\alpha_{\mathrm{s}}\xi_{y}}{2a\sqrt{\rho}}&\frac{\beta_{z}\alpha_{\mathrm{s}}{h_{\xi}}^{2}-\beta_{\xi}\alpha_{\mathrm{s}}\xi_{z}}{2a\sqrt{\rho}}&\frac{\alpha_{\mathrm{f}}{h_{\xi}}^{2}}{2\rho
a^{2}}\bigr{]}.\end{array}$

Details for simplifying the MHD eigenstructure in the curvilinear coordinate
system is given in Appendix B. Similar sets of eigenvectors and eigenvalues
can be defined for $\mathbf{A}_{\eta}$ ($\mathbf{A}_{\zeta}$), by substituting
$\xi_{x}$, $\xi_{y}$, and $\xi_{z}$ with $\eta_{x}$, $\eta_{y}$, and
$\eta_{z}$ ($\zeta_{x}$, $\zeta_{y}$, and $\zeta_{z}$), respectively. Matrices
$\mathbf{P}^{-1}$ and $\mathbf{P}$ for the curvilinear coordinate case are
given as

$\mathbf{P}^{-1}=\begin{bmatrix}J&0&0&0&0&0&0&0\\\
-\frac{u}{\rho}J&\frac{1}{\rho}J&0&0&0&0&0&0\\\
-\frac{v}{\rho}J&0&\frac{1}{\rho}J&0&0&0&0&0\\\
-\frac{w}{\rho}J&0&0&\frac{1}{\rho}J&0&0&0&0\\\ 0&0&0&0&J&0&0&0\\\
0&0&0&0&0&J&0&0\\\ 0&0&0&0&0&0&J&0\\\
\left(\frac{u^{2}+v^{2}+w^{2}}{2}\right)KJ&-uKJ&-vKJ&-wKJ&-B_{x}KJ&-B_{y}KJ&-B_{z}KJ&KJ\end{bmatrix},$
(38)

and

$\mathbf{P}=\begin{bmatrix}\frac{1}{J}&0&0&0&0&0&0&0\\\
\frac{u}{J}&\frac{\rho}{J}&0&0&0&0&0&0\\\
\frac{v}{J}&0&\frac{\rho}{J}&0&0&0&0&0\\\
\frac{w}{J}&0&0&\frac{\rho}{J}&0&0&0&0\\\ 0&0&0&0&\frac{1}{J}&0&0&0\\\
0&0&0&0&0&\frac{1}{J}&0&0\\\ 0&0&0&0&0&0&\frac{1}{J}&0\\\
\left(\frac{u^{2}+v^{2}+w^{2}}{2J}\right)&\frac{\rho u}{J}&\frac{\rho
v}{J}&\frac{\rho
w}{J}&\frac{B_{x}}{J}&\frac{B_{y}}{J}&\frac{B_{z}}{J}&\frac{1}{KJ}\end{bmatrix}$
(39)

where $K=(\gamma-1)$, and $\gamma$ is the specific heat ratio. Using the
derived MHD eigenstructure in the curvilinear coordinates, we can employ the
same procedure as the previous section to apply the characteristic boundary
scheme.

###  4.3 Implementing boundary conditions

Characteristic equations at the boundaries contain eigenvalues of both signs,
showing that characteristic waves are propagating into or out of the
computational domain. Characteristic equations derived in the MHD
characteristic boundary scheme allow us to study each wave separately. The
outgoing waves only depend on the information received from inside the
computational domain and at the boundaries. On the other hand, the incoming
waves carry the information from outside the computational domain. Ignoring
information from outside the computational domain leads to uncertainties in
the approximation of the incoming waves. The number of imposed boundary
conditions should be equal to the number of waves entering the computational
domain. Consider the outlet boundary condition at the right boundary.
Depending on the magnitude of the characteristic eigenvalues, waves
corresponding to the operators $\mathcal{L}_{1}$ to $\mathcal{L}_{3}$ can be
incoming or outgoing. Figure 2 shows the characteristic waves at each
boundary. For super magneto-slow flow, i.e., $u>c_{\mathrm{s}}$, two boundary
conditions can be imposed. On the other hand, only one boundary condition can
be imposed, provided the flow is super Alfvénic, i.e., $u>c_{\mathrm{a}}$.
Finally, no boundary condition can be imposed if the flow is super magneto-
fast, i.e., $u>c_{\mathrm{f}}$, because no wave enters the domain [6].

Fig. 2: Characteristic waves cross the boundaries in the MHD problem.

The time evolution of the primitive variables as a function of
$\mathcal{L}_{i}$ operators for the one-dimensional MHD equations in the
cartesian coordinates can be derived as follows

$\frac{\partial\rho}{\partial
t}+\alpha_{\mathrm{f}}\rho\left(\mathcal{L}_{1}+\mathcal{L}_{8}\right)+\alpha_{\mathrm{s}}\rho\left(\mathcal{L}_{3}+\mathcal{L}_{6}\right)+\mathcal{L}_{4}=0$
(40a) $\frac{\partial u}{\partial
t}+\alpha_{\mathrm{f}}c_{\mathrm{f}}\left(\mathcal{L}_{8}-\mathcal{L}_{1}\right)+\alpha_{\mathrm{s}}c_{\mathrm{s}}\left(\mathcal{L}_{6}-\mathcal{L}_{3}\right)=0,$
(40b) $\frac{\partial v}{\partial
t}+\alpha_{\mathrm{f}}c_{\mathrm{f}}\beta_{y}S\left(\mathcal{L}_{6}-\mathcal{L}_{3}\right)+\alpha_{\mathrm{s}}c_{\mathrm{s}}\beta_{y}S\left(\mathcal{L}_{1}-\mathcal{L}_{8}\right)+\beta_{z}\left(\mathcal{L}_{7}-\mathcal{L}_{2}\right)=0,$
(40c) $\frac{\partial w}{\partial
t}+\alpha_{\mathrm{f}}c_{\mathrm{f}}\beta_{z}S\left(\mathcal{L}_{6}-\mathcal{L}_{3}\right)+\alpha_{\mathrm{s}}c_{\mathrm{s}}\beta_{z}S\left(\mathcal{L}_{1}-\mathcal{L}_{8}\right)+\beta_{y}\left(\mathcal{L}_{2}-\mathcal{L}_{7}\right)=0,$
(40d) $\frac{\partial B_{x}}{\partial t}+\mathcal{L}_{5}=0,$ (40e)
$\frac{\partial B_{y}}{\partial
t}+\alpha_{\mathrm{f}}\sqrt{\rho}\,\beta_{y}\,a\left(-\mathcal{L}_{6}-\mathcal{L}_{3}\right)+\alpha_{\mathrm{s}}\sqrt{\rho}\,\beta_{y}\,a\left(\mathcal{L}_{1}+\mathcal{L}_{8}\right)-\sqrt{\rho}\,\beta_{z}\,S\left(\mathcal{L}_{2}+\mathcal{L}_{7}\right)=0,$
(40f) $\frac{\partial B_{z}}{\partial
t}+\alpha_{\mathrm{f}}\sqrt{\rho}\,\beta_{z}\,a\left(-\mathcal{L}_{6}-\mathcal{L}_{3}\right)+\alpha_{\mathrm{s}}\sqrt{\rho}\,\beta_{z}\,a\left(\mathcal{L}_{1}+\mathcal{L}_{8}\right)+\sqrt{\rho}\,\beta_{y}\,S\left(\mathcal{L}_{2}+\mathcal{L}_{7}\right)=0,$
(40g) and, $\frac{\partial p}{\partial
t}+\alpha_{\mathrm{f}}\,\rho\,a^{2}\left(\mathcal{L}_{1}+\mathcal{L}_{8}\right)+\alpha_{\mathrm{s}}\,\rho\,a^{2}\left(\mathcal{L}_{3}+\mathcal{L}_{6}\right)=0.$
(40h)

Below, the procedure of implementing the non-reflecting inlet/outlet, solid
wall, and constant pressure boundary conditions is discussed.

####  4.3.1 Non-reflecting boundary condition

Applying non-reflective boundary conditions involves setting the magnitude of
the incoming waves to zero [24, 29, 13]. This guarantees that no wave can
enter the computational domain. For example, consider the case of an inlet
boundary condition at the left boundary. The waves associated with the
operators $\mathcal{L}_{4}$ to $\mathcal{L}_{8}$ enter the domain (see Fig.
2), and their amplitude should be set to zero in order to have a non-
reflecting boundary condition. Other waves related to operators
$\mathcal{L}_{1}$ to $\mathcal{L}_{3}$ may be incoming or outgoing waves
depending on their corresponding calculated eigenvalue signs. Non-reflecting
outlet boundary conditions can also be implemented using a similar approach.

####  4.3.2 Slip wall boundary condition

For a slip wall boundary condition, the normal velocity should be set to zero
at the initial condition and remain zero during the simulation. This means
that the condition $\partial u/\partial t=0$ should be met at the boundary. As
Eq. (40b) suggests, for imposing this condition, we set $\mathcal{L}_{6}$
equal to $\mathcal{L}_{3}$ and $\mathcal{L}_{8}$ equal to $\mathcal{L}_{1}$ or
vice versa, depending upon which one leaves or enters the computational
domain. This approach guarantees that the normal velocity component remains
zero and the wall boundary is satisfied.

####  4.3.3 Constant pressure outlet boundary condition

This boundary condition implies that the pressure derivative with respect to
time at the boundary should be zero throughout the simulation. Eq. (40h)
suggests that the time derivative of pressure is controlled by
$\mathcal{L}_{1}$, $\mathcal{L}_{3}$, $\mathcal{L}_{6}$, and
$\mathcal{L}_{8}$. Suppose we have a constant pressure outlet at the right
boundary. In this case, the characteristic waves corresponding to
$\mathcal{L}_{8}$ and $\mathcal{L}_{6}$ are outgoing. Therefore, by setting
$\mathcal{L}_{1}=-\mathcal{L}_{8}$ and $\mathcal{L}_{3}=-\mathcal{L}_{6}$, the
time derivative of the pressure remains zero at the boundary node.

Using the same methodology as the above-mentioned examples, any other kind of
boundary condition can be imposed. The same approach is valid for implementing
different boundary conditions in the curvilinear coordinates.

##  5 Numerical examples

In this section, we first examine the derived characteristic boundary scheme
for the non-magnetic test case, Sod’s shock tube, to verify whether the scheme
works properly for the gasdynamic case. Next, the newly designed test cases
are investigated in order to study the implementation and correctness of the
proposed MHD characteristic boundary scheme for the curvilinear computational
domains.

###  5.1 Sod shock tube test

In this test case, the MHD numerical solver is applied to Sod’s shock tube
problem [28] to assess if the MHD solver can correctly simulate a gasdynamic
case and whether the MHD characteristic boundary scheme can be automatically
reduced to the non-magnetic case. The initial conditions for the left and
right regions are defined in the non-dimensional form as

$\left(\rho,u,p\right)=\begin{cases}(1,0,1)\text{\ \ \ \ \ \ \ \ \ \ \ \ \ \
for $x<0.5$},\\\ (0.125,0,0.1)\ \ \text{\ \ for $x>0.5$},\end{cases}$ (41)

with the heat capacity ratio $\gamma=1.4$. The fifth-order WENO and fourth-
order Runge-Kutta schemes are used for the spatial discretization and time
marching, respectively. The CFL number is $0.5$ and a solid wall boundary
condition is applied at both ends. The computational domain consists of $600$
non-uniform stretched grid points. The size of the computational domain is
$\xi\in[0,1]$, and the grid coordinates in the physical domain are given by

$x=\left(x_{\mathrm{max}}-x_{\mathrm{min}}\right)\frac{\exp{\left(\delta_{x}\,\xi\right)}-1}{\exp\left(\delta_{x}\right)-1}+x_{\mathrm{min}},$
(42)

where $\delta_{x}$ is the stretch factor in the $x$-direction, which is set to
$\delta_{x}=0.8$ for this test case. Parameters $x_{\mathrm{min}}$ and
$x_{\mathrm{max}}$ denote the minimum and maximum values of the physical
domain, respectively, with values $x_{\mathrm{min}}=0$ and
$x_{\mathrm{max}}=1$. Figure 3 shows the distribution of the physical domain
points against the computational ones. It can be observed that the points are
more compressed near $x=0$ and become more stretched as approaching $x=1$.

Fig. 3: Distribution of the physical domain points, $x$, against the
computational ones, $\xi$, in the non-uniform stretched computational domain,
using $20$ grid points with $\delta_{x}=0.8$.

This problem is simulated using two numerical solvers, the implemented MHD
solver as well as the verified in-house Euler solver with the characteristic
boundary scheme. In the MHD solver, all the magnetic components are set to
zero, and the density, velocity, and pressure fields are initialized according
to the initial conditions of the Sod shock tube problem. Setting the magnetic
field to zero results in the disappearance of the slow magnetosonic and Alfvén
waves. Furthermore, the fast magnetosonic waves are reduced to acoustic ones.

Figure 5.1 shows the time evolution of the density and pressure in the
interval $t\in[0,3]$. As can be observed, results from the Euler and MHD
solvers are in excellent agreement. Moreover, propagating waves in the
computational domain have the same patterns and the waves are correctly
reflected from the wall. This observation suggests that the MHD characteristic
boundary scheme works properly for the $\mathbf{B}=\textbf{0}$ case and the
eigensystem is well defined for the non-magnetic case.

(a)  (c)  (d)  Fig. 4: Density and pressure contours using Euler and MHD
solvers on a stretched mesh with $600$ grid points and $\text{CFL}=0.5$ for
$t=0$ to $3$. Wall boundary condition is applied at both ends. (a) and (b)
show the density and pressure contours, respectively, using the Euler solver.
(c) and (d) depict the density and pressure contours, respectively, using the
MHD solver.

For a better comparison, the pressure and density fields obtained from the
Euler and the MHD solvers at $t=0.9375$ 111This time is chosen for the sake of
a better presentation of the results. are shown in the same plot in Fig. 5.1.
The results from the MHD solver are in close agreement with the Euler solver,
which shows that the implemented MHD characteristic boundary scheme can
correctly simulate the characteristic waves interaction with the wall
boundary.

(a)  Fig. 5: Density and pressure results using Euler and MHD solvers on a
stretched mesh with $600$ grid points and $\text{CFL}=0.5$ at $t=0.9375$. (a)
and (b) show the density and pressure results, respectively. —, using Euler
solver; $\bullet$, using MHD solver. For the MHD solver, every 6th point of
the results is shown.

###  5.2 Smooth Alfvén test with non-reflecting boundary condition

The smooth Alfvén test case is widely used in the literature to determine the
order of accuracy of numerical schemes for the MHD problems with smooth flows.
The initial condition is given as

$\left(\rho,u,v,w,p,B_{x},B_{y},B_{z}\right)=\left(1,0,0.1\sin\left(2\pi
x\right),0.1,0.1\cos\left(2\pi x\right),1,0.1\sin\left(2\pi
x\right),0.1\cos\left(2\pi x\right)\right).$ (43)

In previous studies, periodic boundary conditions were applied along both
$x$\- and $y$-directions. Here, we employ the non-reflecting boundary
condition along the $x$-direction and the periodic one along the
$y$-direction. The analytical solution of the smooth Alfvén problem includes
the Alfvén wave moving along the $x$-direction with a speed of $1$ and with
constant density during the simulation. The computational and physical domains
extend to $[0,1]\times[0,1]$; therefore by imposing a non-reflecting boundary
condition along the $x$-direction after $t=1$, the Alfvén wave leaves the
computational domain. As a result, the initial non-uniformities leave the
domain and the magnetic field becomes uniform with
$\mathbf{B}=\left(1,0,0.1\right)$. Hence, the final solution for $t\geq 1$ is
given as

$\left(\rho,u,v,w,p,B_{x},B_{y},B_{z}\right)=\left(1,0,0,0.1,0.1,1,0,0.1\right).$
(44)

In this test case, three computational domains with $128\times 128$ grid size
are selected to investigate the efficiency of the non-reflecting
characteristic boundary scheme. The first computational domain includes
uniform grid points in both $x$\- and $y$-directions, whereas the second one
consists of non-uniform stretched grid points along the $x$-direction using
Eq. (42) and uniform grid points in the $y$-direction. The third computational
domain is set to be $[\xi,\eta]^{2}\in[0,1]\times[0,1]$ with the perturbed
$x$-grid lines according to the following mapping function [5]

$x=\xi+\alpha_{x}\sin\left(2\pi\,\eta\,\beta_{x}\right),$ (45)

where $\alpha_{x}$ shows the magnitude of the perturbation and $\beta_{x}$ is
the wave number of the perturbation [5] with uniform $y$-grid points. In this
study, parameters $\alpha_{x}$ and $\beta_{x}$ are defined the same as the
ones mentioned by Christlieb et al. [5], where $\alpha_{x}$ and $\beta_{x}$
are equal to $0.01$ and $2$, respectively. These three computational domains
are shown in Fig. 5.2.

As explained in section 4.3.1, a non-reflecting boundary condition is imposed
by setting the amplitude of the incoming waves to zero [29]. In this test
case, we investigate the implementation and accuracy of the non-reflecting
boundary condition using the MHD characteristic boundary scheme. The sixth-
order compact and fourth-order Runge-Kutta schemes are used for the spatial
discretization and time marching, respectively, with CFL$=0.6$.

(a)  (c)  Fig. 6: Three computational domains for solving the two-dimensional
smooth Alfvén test with non-reflecting boundaries along $x$-direction. (a)
shows uniform grid points in both $x$\- and $y$-directions. Computational
domain (b) has a non-uniform stretched grid points in $x$-direction using Eq.
(42), and uniform grid points in $y$-direction. In computational domain (c),
grid points along $x$-direction are perturbed according to Eq. (45) while $y$
grid points remain uniform. Every 4th node is shown.

The calculated Root Mean Square (RMS) of the error for the magnetic field
components is less than $0.05\%$ for the uniform computational domain. The RMS
of the errors is calculated using

$\text{RMSE}=\sqrt{\frac{\Sigma_{i=1}^{N}\left(f_{\mathrm{anal}_{i}}-f_{\mathrm{num}_{i}}\right)^{2}}{N}},$
(46)

where $f_{\mathrm{anal}_{i}}$ and $f_{\mathrm{num}_{i}}$ show the analytical
and numerical results at grid point $i$, respectively, and $N$ is the total
number of grid points. For the stretched grid, the results perfectly followed
the analytical ones and the error is almost $0.1\%$ for the magnetic field
components. Small discrepancies are observed compared to the analytical
solution for the curvilinear computational domain, since the characteristic
boundary scheme is derived while assuming the flow to be locally one-
dimensional. This assumption is not valid for the curvilinear mesh and thus,
these errors are inevitable. The calculated RMSE in the curvilinear domain for
$B_{x}$, $B_{y}$, $B_{z}$, and $p$ are $0.23\%$, $0.21\%$, $0.22\%$, and
$0.27\%$, respectively. Calculated errors for the three computational domains
are reported in Table (1). Furthermore, it is also beneficial to calculate the
$L^{\infty}$ error for the curvilinear computational domain since it shows the
maximum error and evaluate the scheme’s performance on the grid points with
the greatest distortion. The obtained $L^{\infty}$ error for $B_{x}$, $B_{y}$,
$B_{z}$, and $p$ on the curvilinear domain are 0.44%, 0.47%, 0.33%, and 0.29%,
respectively. As expected, the calculated $L^{\infty}$ errors are higher than
$L^{2}$ ones. However, since the computed errors are under 0.5%, we can
conclude that the scheme works well on boundary grids with the greatest
misalignment as well.

Table 1:  Calculated RMSE for the two-dimensional smooth Alfvén wave test case with the non-reflecting boundaries along $x$-direction, using three computational domains. Computational domain | Error in $p$ | Error in $B_{x}$ | Error in $B_{y}$ | Error in $B_{z}$
---|---|---|---|---
(a) | $\leq 0.06\%$ | $\leq 0.05\%$ | $\leq 0.04\%$ | $\leq 0.05\%$
(b) | $0.13\%$ | $0.1\%$ | $0.11\%$ | $0.1\%$
(c) | $0.27\%$ | $0.23\%$ | $0.21\%$ | $0.22\%$

Figure 5.2 shows the comparison of the magnetic field and pressure results
using three computational domains. Results from the uniform and non-uniform
computational domains are indistinguishable on the same plot, while for the
curvilinear case, small discrepancies are observed.

(a)  (c)  (d)  Fig. 7: Comparing the magnetic field and pressure results using
three computational domains for the two-dimensional smooth Alfvén problem with
non-reflecting boundaries along $x$ and periodic boundary along $y$, at $t=1$
and $y=0.5$. —, uniform computational domain; $\times$, non-uniform stretched
computational domain; $\bullet$, curvilinear computational domain. For the
stretch and the curvilinear grids, every 5th and 4th point of the results is
shown, respectively.

###  5.3 Gaussian pulse test with wall boundary condition

This one-dimensional test case is designed to investigate the wall boundary
condition for the MHD problem using the implemented characteristic boundary
scheme. The initial condition reads as

$(\rho,u,v,B_{x},B_{y},p)=\left(1.0,0,0,0.75,1+\epsilon\exp\left[\frac{-x^{2}}{b^{2}}\right],1.1\left(1+\epsilon\exp\left[\frac{-x^{2}}{b^{2}}\right]\right)\right),$
(52)

where $\epsilon=0.1$, and $b=1.2$. As clearly seen in Fig. 5.3, the initial
condition consists of a Gaussian pulse perturbation applied to the pressure
and the transverse magnetic field at $x=0$. The no-penetration wall boundary
condition is applied at both ends.

(a)  (c)  (d)  Fig. 8: The one-dimensional Gaussian pulse problem using $400$
grid points with wall boundaries, at $t=0$. (a) density, (b) pressure, (c)
transverse component of the magnetic field, and (d) velocity.

The normal velocity is forced to zero at the wall. According to Eq. (40b), we
need to set $\mathcal{L}_{6}$ equal to $\mathcal{L}_{3}$ and $\mathcal{L}_{8}$
to $\mathcal{L}_{1}$ at the left boundary and the other way around for the
right boundary to satisfy the wall boundary condition. This procedure ensures
that $\partial u/\partial t$ remains zero during the simulation, which is
essential for this type of boundary condition.

The simulation is conducted in the non-uniform stretched grid according to Eq.
(42) mapping function, with $400$ grid points in the computational domain,
$\xi\in[-20,20]$, and $\text{CFL}=0.6$. The sixth-order compact scheme and the
fourth-order Runge-Kutta scheme are used for the spatial discretization and
time marching, respectively. Moreover, this test case is simulated in the
uniform computational domain to compare the results with its non-uniform
counterpart.

Figure (5.3) show the evolution of the density, pressure, transverse magnetic,
and velocity fields before and after the interaction of the waves with the
walls. As can be observed, there is a good agreement between the uniform and
non-uniform results. This suggests that the characteristic boundary scheme
works properly, independent of the grid configuration.

Since in the initial condition, the flow speed is assumed to be zero, the
results are symmetric with respect to the origin, $x=0$, for the pressure,
density, and the transverse magnetic fields. For the velocity field, the left-
and right-going characteristic waves are observed having the same speed,
albeit with different signs.

Figure 5.3 displays the results at $t=5$ in black, before the interaction of
the waves with the wall boundaries. At $t=5$, the propagation of the left and
right, fast and slow magnetosonic waves can be observed. There is a pressure
bump at the origin at $t=0$ and, hence, the density value should decrease at
$x=0$ to compensate the pressure disturbance at this point. Figure 5.3 shows
the results as the left and right going waves interact with the wall at
$t=11.4$ in red. The magnitude of the pressure increases at both ends due to
the presence of the wall. Figure 5.3 depicts the reflected waves from the wall
boundaries at $t=15$ in blue. The reflected waves are symmetric with respect
to the origin, as expected. Moreover, the reflected density and pressure peaks
have the same values as those before their interaction with the wall. These
observations validate the implemented wall boundary condition.

(a)  (c)  (d)  Fig. 9: The one-dimensional Gaussian pulse problem with wall
boundaries; —, non-uniform computational domain; $\bullet$, uniform
computational domain with $400$ grid points at $t=5$ before waves interact
with the walls (in black), $t=11.4$ as the waves interact with the walls (in
red), and $t=15$ after the waves reflected from the walls (in blue). Every 5th
point of the results is shown for the uniform case. (a) density, (b) pressure,
(c) transverse component of the magnetic field, and (d) velocity.

###  5.4 Magnetosonic wave propagation test

####  5.4.1 One-dimensional case

The one-dimensional magnetosonic wave propagation test case is designed to
investigate the single wave inlet boundary condition. This test case
demonstrates the procedure of using the characteristic boundary scheme to
inject a single magnetosonic wave into the domain at a boundary. Subsequently,
the evolution of the injected characteristic wave into the computational
domain can be studied. The initial condition, whose uniformity helps us to
better verify the time evolution of the injected characteristic wave, is given
as

$(\rho,u,v,B_{x},B_{y},p)=\left(1.0,0,0,0.75,1.0,0.5\right),$ (53)

with the single fast magnetosonic wave inlet and non-reflecting outlet
boundary conditions applied to the left and right boundaries, respectively.
For this test case, the specific heat value is assumed to be $\gamma=1.4$.

In order to have a single fast characteristic wave entering the domain, we
should set the magnitude of all the incoming waves to zero, except the fast
one. Here, we assume a Gaussian pulse is moving towards the left boundary with
the fast magnetosonic speed. To implement the fast magnetosonic inlet boundary
condition, we defined the $\mathcal{L}_{8}$ operator, the characteristic
operator with respect to the fast characteristic wave, as:

$\mathcal{L}_{8}=-\epsilon\exp\left[\frac{-\left(c_{\mathrm{f}}\,t-1\right)^{2}}{2b^{2}}\right],$
(54)

where variables $\epsilon$ and $b$ are assumed to be 0.5 and 0.2,
respectively. Furthermore, $c_{\mathrm{f}}$ shows the velocity of the injected
fast wave. Since $\mathcal{L}_{8}$ is the only incoming wave, the time
evolution of the pressure can be written as

$\frac{\partial p}{\partial
t}+\alpha_{\mathrm{f}}\,\rho\,a^{2}\mathcal{L}_{8}=0.$ (55)

Therefore, evolution of the pressure is given by the error function, i.e., the
integral of the Gaussian pulse function. Similarly, the time evolution of the
other primitive variables are

$\frac{\partial\rho}{\partial t}+\alpha_{\mathrm{f}}\,\rho\mathcal{L}_{8}=0,$
(56) $\frac{\partial u}{\partial
t}+\alpha_{\mathrm{f}}\,c_{\mathrm{f}}\mathcal{L}_{8}=0,$ (57)

and,

$\frac{\partial B_{y}}{\partial
t}+\alpha_{\mathrm{s}}\,c_{\mathrm{s}}\,\sqrt{\rho}\,\beta_{y}\,a\mathcal{L}_{8}=0.$
(58)

According to Eqs. (56) to (58), we expect the results for the density,
velocity, and the transverse magnetic fields to be in the form of the error
function as well.

This test is also repeated for the case of injecting a Gaussian pulse with the
slow magnetosonic speed into the computational domain at the left boundary,
with the same initial condition. Analogous to the previous case, we again set
the magnitude of all the incoming waves to zero, with the exception of the
slow one. Hence, the $\mathcal{L}_{6}$ operator is defined as

$\mathcal{L}_{6}=-\epsilon\exp\left[\frac{-\left(c_{\mathrm{s}}\,t-1\right)^{2}}{2b^{2}}\right],$
(59)

where variables $\epsilon$ and $b$ are assumed to be 0.5 and 0.2,
respectively.

The simulation of both the slow and fast characteristic wave injections are
conducted on a computational domain with the non-uniform stretched grid points
constructed according to Eq. (42), with $200$ grid points and
$\text{CFL}=0.5$, for $x\in[0,5]$. Results for these two cases are shown in
Fig. 5.4.1.

(a)  (c)  (d)  Fig. 10: The magnetosonic waves propagation problem on a non-
uniform stretched grid with $200$ points at $t=3.5$. (a) density, (b)
pressure, (c) transverse component of the magnetic field, and (d) velocity. —,
fast magnetosonic wave propagation; $\bullet$, slow magnetosonic wave
propagation.

As expected, for the fast case, the injected disturbance travels faster in the
domain, and therefore the initial density, velocity, pressure, and magnetic
fields evolve more rapidly. The calculated mean speed value for the slow and
fast magnetosonic waves are $0.435$ and $1.439$, respectively. Hence, at
$t=3.5$, we expect the position lag between the onset of the error functions
to be approximately $3.5$, which can be observed in Fig. 5.4.1 as well.

The difference exhibited in the amplitudes of the results are due to the
difference in the dependence of the time derivatives of the primitive
variables on the amplitudes of the waves. Furthermore, as already seen in Eq.
(40f), the coefficients of the $\mathcal{L}_{6}$ and $\mathcal{L}_{8}$
operators have different signs. Therefore, the signs of the induced error
functions for the transverse magnetic fields are opposite.

The accuracy of the results can also be verified by studying the density and
magnetic fields reaction to the perturbation injected into the domain. As
mentioned earlier, the density and magnetic fields are in phase for the fast
magnetosonic wave. This behaviour can be observed in Fig. 5.4.1 as they both
have the same trend. However, for the slow magnetosonic wave, density and
magnetic fields are out of phase. This is why their behaviour is opposite for
the slow magnetosonic test case.

We also studied this test case by injecting an Alfvén wave at the left
boundary. As the Alfvén wave is incompressible, we expect its effect on
density to be null. In order to apply this boundary condition, we define the
$\mathcal{L}_{7}$ operator as

$\mathcal{L}_{7}=-\epsilon\exp\left[\frac{-\left(c_{\mathrm{a}}\,t-1\right)^{2}}{2b^{2}}\right],$
(60)

where variables $\epsilon$ and $b$ are set similarly as the previous test.
Figure 5.4.1 shows the results of this test case. In the time evolution
relations for the density, velocity, and pressure, Eqs. (40a), (40b), and
(40h), the $\mathcal{L}_{7}$ operator does not appear, and, hence, has no
effect on these variables. On the other hand, the time evolution of the
transverse magnetic field can be recast as

$\frac{\partial B_{y}}{\partial
t}-\sqrt{\rho}\,\beta_{z}\,S\left(\mathcal{L}_{2}+\mathcal{L}_{7}\right)=0.$
(61)

Therefore, the $\mathcal{L}_{7}$ operator affects the transverse magnetic
field but since $B_{z}=0$, the coefficient of the Alfvén characteristic waves
becomes zero. Consequently, as shown in Fig. 5.4.1, the transverse magnetic
field does not change at all.

(a)  (c)  (d)  Fig. 11: The Alfvén wave propagation problem on a non-uniform
stretched grid with $200$ points at $t=1.0$. (a) density, (b) pressure, (c)
transverse component of the magnetic field, and (d) velocity.

####  5.4.2 Two-dimensional case

The fast magnetosonic wave propagation test case is also conducted on a two-
dimensional computational domain similar to Fig. 5.2c with
$[x,y]\in[0,5]\times[0,5]$ and the same initial condition as Eq. (53). At the
left boundary along the $x$-direction, a fast magnetosonic wave is injected
into the domain, while the right boundary is set to be a non-reflecting
outlet. Periodic boundaries are applied along the $y$-direction. The
simulation is performed on $200\times 100$ grid points and the CFL number is
$0.6$.

Results from the perturbed two-dimensional computational domain at $y=2.5$ are
compared with those from the one-dimensional case, shown in Fig. 5.4.2 at
$t=1.86$ 222This time is chosen for the sake of a better presentation of the
results.. The calculated RMSE for the density, pressure, velocity, and
transverse magnetic fields from the two-dimensional simulation relative to the
one-dimensional case are $0.30\%$, $0.21\%$, $0.42\%$, and $0.41\%$,
respectively. Therefore, we conclude that the single wave injection boundary
condition works properly on the curvilinear domain as well. The primary
sources of the captured errors for this test case are the perturbed boundary
points along the $x$-direction, which are functions of $\xi$ and $\eta$. This
computational domain is not consistent with the locally one-dimensional
assumption that we made while implementing the characteristic boundary scheme.
Therefore, these errors are prone to occur.

(a)  (c)  (d)  Fig. 12: The two-dimensional fast magnetosonic wave propagation
problem using $200\times 100$ grid points on a perturbed computational domain
at $t=1.86$. (a) density, (b) pressure, (c) transverse component of the
magnetic field, and (d) velocity. —, results from two-dimensional test case at
$y=2.5$; $\bullet$, results from one-dimensional case. For the one-dimensional
test case, every 4th point is shown.

###  5.5 Brio–Wu shock tube test with non-reflecting boundary condition

This test case was first proposed by Brio and Wu [3]. The initial condition is
given as:

$\left(\rho,u,v,w,p,B_{x},B_{y},B_{z}\right)=\begin{cases}(1,0,0,0,1,0.75,+1,0)\
\ \ \ \ \ \ \ \ \ \ \ \text{\ \ for $x<0$},\\\ (0.125,0,0,0,0.1,0.75,-1,0)\ \
\ \ \text{for $x>0$},\end{cases}$ (62)

with the heat capacity ratio $\gamma=2$. The Brio–Wu problem is a standard
test for MHD codes since it accurately represents the shocks, rarefactions,
contact discontinuities, and compound structure of the MHD problem. This test
case exhibits strong shocks and discontinuities in the presence of a magnetic
field, and thus comprises a rigorous test case for the proposed characteristic
boundary scheme for MHD problems that are of interest to applications such as
fusion energy.

Due to the existence of sharp discontinuities and shocks, it is necessary to
apply Total Variation Diminishing (TVD) schemes for temporal integration and
spatial discretization. That being the case, in this problem, the third-order
TVD Runge-Kutta [12] and fifth-order WENO methods have been utilized for the
time integral and spatial discretization, respectively. Furthermore, non-
reflecting boundaries are applied at the right and left boundaries for this
case. The non-reflecting boundary condition will help us assess whether the
implemented characteristic boundary method can efficiently distinguish the
incoming and outgoing waves at the inlet and outlet, especially when shocks
and jump discontinuities reach the boundaries. Furthermore, the validity of
the derived eigenvectors and eigenvalues can be evaluated for problems with
discontinuities. The test case is simulated on the computational domain
consisting of $800$ non-uniform stretched grid points according to Eq. (42)
with $x\in[-1,1]$. The CFL number is set to $0.6$ with the simulation time
$t=0$ to $t=3$.

The solution of the Brio–Wu shock tube consists of a right-traveling magneto-
fast rarefaction wave, which is followed by a magneto-slow shock wave and a
contact discontinuity. Moreover, two waves, i.e., magneto-fast rarefaction and
magneto-slow compound wave, propagate to the left. The compound wave contains
an intermediate shock followed by a slow rarefaction wave. Figure (13) depicts
the evolution of the velocity field at six different time steps, $t=0$,
$0.121$, $0.229$, $0.602$, $0.940$, and $3$. At $t=0.121$, both right- and
left-traveling shock and rarefaction waves can be observed. At $t=0.229$, the
fast rarefaction wave propagating to the right gets closer to the boundary.
According to Fig. (13), at $t=0.602$, the right-traveling rarefaction wave has
already left the right boundary without any reflection. Moreover, the other
rarefaction wave moving towards the left boundary is also leaving the domain.
At $t=0.940$, the right-traveling slow shock wave has completely left the
computational domain and the left-traveling rarefaction wave is also on the
verge of crossing the left boundary entirely. At the last time step, $t=3$,
the left-traveling compound wave has gotten closer to the left boundary, as we
expected, and the left-traveling rarefaction wave has completely exited the
domain. Figure (14) represents the evolution of the transverse magnetic field
at the same time steps. Similarly, it can be observed that the right- and
left-traveling rarefactions and shocks leave the domain at the left and right
non-reflecting boundaries without causing any numerical instability due to the
interaction of discontinuities with the boundaries. Furthermore, it can be
concluded that the proposed characteristic boundary scheme performs robustly
in the presence of shocks as well, and the derived eigenstructure remains
valid.

Fig. 13: The velocity field of the Brio–Wu shock tube test case with non-
reflecting boundaries on a stretched mesh with $800$ grid points and
$\text{CFL}=0.6$ at six different time steps, $t=0$, $0.121$, $0.229$,
$0.602$, $0.940$, and $3$. Fig. 14: The transverse magnetic field of the
Brio–Wu shock tube test case with non-reflecting boundaries on a stretched
mesh with $800$ grid points and $\text{CFL}=0.6$ at six different time steps,
$t=0$, $0.121$, $0.229$, $0.602$, $0.940$, and $3$.

##  6 Conclusion

Imposing proper boundary conditions is of crucial importance for the accurate
simulation of plasma dynamics. In this study, we introduce characteristic
boundary conditions for the numerical solution of the MHD equations on multi-
dimensional, curvilinear grids. A parallel, high-order code, with sixth-order
accuracy in space and fourth-order accuracy in time, was developed for the
simulation of the compressible MHD equations. A fifth-order WENO scheme was
also implemented to capture discontinuities such as shock waves.

The correctness of the derived scheme was first verified for the Sod’s shock
tube problem, a non-magnetic test case. Then, a set of test cases were
designed to investigate the accuracy of the proposed scheme, in particular
non-reflecting boundary conditions which are of crucial importance for the
accurate simulations of instabilities. We have also shown how to inject
single-mode waves into the computational domain. Such a boundary condition can
be used to independently study the dynamics of each wave mode. In conclusion,
the studied test cases showed the promise and correctness of the implemented
MHD characteristic boundary scheme for one- and two-dimensional problems in
both cartesian and curvilinear coordinate systems.

## Acknowledgments

The authors thank Andrew Higgins for his invaluable suggestions and comments
to improve this manuscript. This study was supported by the Natural Science
and Engineering Research Council of Canada (NSERC) through an NSERC Discovery
Grant.

## Appendix A. Convergence study

## Two-Dimensional Smooth Alfvén Test

This problem is widely used to verify the accuracy order of the numerical
scheme for MHD problems with smooth flows. The capability of Alfvén waves to
move over long times and distances in numerical solvers is important
especially in turbulence simulations [9]. Alfvén waves can be damped as a
result of numerical dissipation in the scheme. If the code fails to capture
these waves, the turbulence behaviour cannot be correctly studied, as MHD
turbulence is mainly sustained by Alfvén waves[1]. The initial condition is
given as:

$\left(\rho,u,v,w,p,B_{x},B_{y},B_{z}\right)=\left(1,0,0.1\sin\left(2\pi
x\right),0.1,0.1\cos\left(2\pi x\right),1,0.1\sin\left(2\pi
x\right),0.1\cos\left(2\pi x\right)\right),$ (63)

with the periodic boundary condition along the $x$\- and $y$-directions. The
computational domain is set to be $[\xi,\eta]^{2}\in[0,1]\times[0,1]$ with the
perturbed grid lines according to the mapping functions [5]

$x=\xi+\alpha_{x}\sin\left(2\pi\,\eta\,\beta_{x}\right),$ (64)

and

$y=\eta+\alpha_{y}\sin\left(2\pi\,\xi\,\beta_{y}\right),$ (65)

where $\alpha_{x}$ and $\alpha_{y}$ show the magnitude of the perturbation and
$\beta_{x}$ and $\beta_{y}$ are the wave numbers of the perturbation [5]. In
this study, these parameters are taken to be the same as the ones mentioned by
Christlieb et al. [5], where $\alpha_{x}$, $\alpha_{y}$, $\beta_{x}$, and
$\beta_{y}$ are equal to $0.01$, $0.02$, $2$, and $4$, respectively. The
computational domain is shown in Fig. (15).

Fig. 15: Computational domain for the two-dimensional smooth Alfvén wave
problem with periodic boundary conditions along the $x$\- and $y$-directions.
The grid size of $32\times 32$ is shown with
$[\xi,\eta]^{2}\in[0,1]\times[0,1]$.

The exact solution is an Alfvén wave propagating along the $x$-direction with
a speed of 1. Therefore, the initial circularly polarized Alfvén wave
propagates across a periodic domain. The numerical results, using both spatial
discretization schemes, the central compact and the WENO schemes, are compared
to the analytical solution at $t=5$. Different grid resolutions with the
constant $\text{CFL}=0.5$ were used for calculating the order of accuracy, as
reported in Tables (2) and (3). In Tables (2) and (3), the $L^{\infty}$ error
values are shown as well as the estimated convergence rates. For each
resolution, the order of accuracy was calculated with respect to the previous
courser one. As can be concluded from Tables 2 and 3, fifth-order WENO and
sixth-order compact schemes are almost fourth-order and fifth-order globally
accurate, respectively.

Table 2: $L^{\infty}$ errors for two-dimensional smooth Alfvén wave problem
using the fifth-order WENO scheme. WENO scheme
---
Mesh |  Error in $\mathbf{B}$ |  Order
32 $\times$ 32 |  5.4128-3 |  ——
64 $\times$ 64 |  4.3418e-4 |  3.64
128 $\times$ 128 |  2.8093e-5 |  3.95
256 $\times$ 256 |  1.7680e-6 |  3.99
Table 3: $L^{\infty}$ errors for two-dimensional smooth Alfvén wave problem
using the sixth-order compact scheme. Compact scheme
---
Mesh |  Error in $\mathbf{B}$ |  Order
32 $\times$ 32 |  4.9879e-3 |  ——
64 $\times$ 64 |  1.9298e-4 |  4.66
128 $\times$ 128 |  6.2869e-6 |  4.94
256 $\times$ 256 |  1.9920e-7 |  4.98

## Appendix B. Simplifying MHD eigenstructure

As was highlighted before, three types of waves termed fast, intermediate
(Alfvén), and slow are formed according to their characteristic speeds in the
MHD system. Depending on the magnitude and direction of the magnetic field,
these waves may coincide with each other. Therefore, the MHD equations may
involve multiplicities of eigenvalues, provided the modes coincide. For
example, this phenomenon may occur in the following conditions for the one-
dimensional case [4, 3]: 1) When $B_{x}=0$ and
$c_{\mathrm{s}}=c_{\mathrm{a}}=0$. Therefore, the eigenvalue $u$ becomes of
multiplicity five. 2) When ${B_{y}}^{2}+{B_{z}}^{2}=0$,
${c_{\mathrm{s}}}^{2}=\min(a^{2},{b_{x}}^{2})$, and
${c_{\mathrm{f}}}^{2}=\max(a^{2},b_{x}^{2})$. Thus for the case
$a^{2}\neq{b_{x}}^{2}$, either if ${c_{\mathrm{s}}}^{2}={b_{x}}^{2}$ or
${c_{\mathrm{f}}}^{2}={b_{x}}^{2}$, the multiplicity of $u\pm c_{\mathrm{a}}$
becomes two. Also, for the case
${c_{\mathrm{s}}}^{2}={c_{\mathrm{f}}}^{2}={b_{x}}^{2}$, the multiplicity of
$u\pm c_{\mathrm{a}}$ changes to three. Consequently, the MHD governing
equations do not form the strictly hyperbolic system of equations. Owing to
this property, we should make sure that the eigensystem of the MHD equations
remains well-defined for all cases.

The first study of the MHD eigenvectors was carried out by Jeffrey and Tanuiti
[15] and was continued by others. The derived set of eigenvectors by Jeffrey
and Tanuiti does not form a complete set near the points $B_{x}=0$ or
${B_{y}}^{2}+{B_{z}}^{2}=0$. Since in the proximity of the abovementioned
conditions, this set of eigenvectors is not well defined as its columns become
singular [4]. To address this issue, a set of normalization parameters was
proposed by Brio and Wu [4], and a new set of eigenvectors was derived for the
cartesian coordinate system. These normalization parameters are
$\alpha_{\mathrm{s}}$, $\alpha_{\mathrm{f}}$, and $\beta_{x(y,z)}$, defined in
Eq. (35) [26]. These normalization parameters were proposed by others [26, 31,
8] as well. For these parameters, the following relations are held

${\alpha_{\mathrm{f}}}^{2}+{\alpha_{\mathrm{s}}}^{2}=1,\ \
{\alpha_{\mathrm{f}}}^{2}{c_{\mathrm{f}}}^{2}+{\alpha_{\mathrm{s}}}^{2}{c_{\mathrm{s}}}^{2}=a^{2},\
\
{\alpha_{\mathrm{s}}}^{2}{c_{\mathrm{f}}}^{2}+{\alpha_{\mathrm{f}}}^{2}{c_{\mathrm{s}}}^{2}=b^{2},\
\
\alpha_{\mathrm{f}}\alpha_{\mathrm{s}}=\frac{a\,b_{\perp}}{{c_{\mathrm{f}}}^{2}-{c_{\mathrm{s}}}^{2}},$
(78)

where
${c_{\mathrm{f},\mathrm{s}}}^{2}=1/2\left[a^{2}+b^{2}\pm\sqrt{\left(a^{2}+b^{2}\right)^{2}-4\,a^{2}\,{b_{x}}^{2}}\right]$,
$a^{2}=\gamma p/\rho$, $b={b_{x}}^{2}+{b_{y}}^{2}+{b_{z}}^{2}$, and
${b_{\perp}}^{2}={b_{y}}^{2}+{b_{z}}^{2}$. Moreover, it can be shown that the
square values of the fast and slow magnetosonic waves (${c_{\mathrm{f}}}^{2}$
and ${c_{\mathrm{s}}}^{2}$) are the bigger and smaller roots of the following
equation [26]:

$c^{4}-(a^{2}+b^{2})c^{2}+a^{2}{b_{x}}^{2}=0,$ (79)

where $c^{2}$ is an independent variable. Using Eq. (79), we can easily prove
that the following relations are satisfied:

${c_{\mathrm{f}}}^{2}+{c_{\mathrm{s}}}^{2}=a^{2}+b^{2},\ \
c_{\mathrm{f}}\,c_{\mathrm{s}}=a\,|b_{x}|.$ (80)

Eqs. (78) and (80) are widely used for simplifying the eigensystem of the MHD
equations in the cartesian coordinate system [26]. In order to extend the MHD
eigenstructure to the general three-dimensional curvilinear coordinate system,
a new set of normalization parameters should be derived. In this study, these
parameters are defined similar to the cartesian coordinates by substituting
variables $a^{2}$, $b^{2}$, $b_{\perp}$, $c_{\mathrm{f}}$, $c_{\mathrm{s}}$,
and $c_{\mathrm{a}}$ with the revised ones derived for the curvilinear case
(see Eqs. (31) to (33)). Equations (78) to (80) are also held for the new
defined normalization parameters in the general curvilinear coordinate system.
In this work, these relations are widely used to simplify the set of
calculated eigenvectors as much as possible and prevent any indeterminacy in
the proposed set of eigensystem.

## References

  * [1] D. S. Balsara.  Multidimensional reimann problem with self-similar internal structure. part 1 - application to hyperbolic conservation laws on structured meshes.  Journal of Computational Physics, 277:163–200, 2014.
  * [2] A. Bayliss and E. Turkel.  Far field boundary conditions for compressible flows.  Journal of Computational Physics, 48(2):182–199, 1982.
  * [3] M. Brio and C. C. Wu.  Characteristic fields for the equations of magnetohydrodynamics.  Contemporary Mathematics, 60:19–23, 1987.
  * [4] M. Brio and C. C. Wu.  An upwind differencing scheme for the equations of ideal magnetohydrodynamics.  Journal of Computational Physics, 75(2):400–422, 1988.
  * [5] A. J. Christlieb, X. Feng, Y. Jiang, and Q. Tang.  A high-order finite difference WENO scheme for ideal magnetohydrodynamics on curvilinear meshes.  SIAM Journal on Scientific Computing, 40(4):A2631–A2666, 2017.
  * [6] A. Cimino, G. Krause, S. Elaskar, and A. Costa.  Characteristic boundary conditions for magnetohydrodynamics: The Brio–Wu shock tube.  Computers and Fluids, 127:194–210, 2016.
  * [7] S. D. Conte and C. Boor.  Elementary Numerical Analysis: An Algorithmic Approach, third ed. McGraw-Hill, Inc, 1980.
  * [8] W. Dai and P. R. Woodward.  Extension of the piecewise-parabolic method to multidimensional ideal magnetohydrodynamics.  Journal of Computational Physics, 115(2):485–514, 1994.
  * [9] D. Derigs, A. R. Winters, G. J. Gassner, and S. Walch.  A novel high-order, entropy stable, 3D AMR MHD solver with guaranteed positive pressure.  Journal of Computational Physics, 317:223–256, 2016.
  * [10] B. Engquist and A. Majda.  Absorbing boundary conditions for the numerical solution of waves.  Mathematics of Computation, 31:629–651, 1977.
  * [11] B. Engquist and A. Majda.  Radiation boundary conditions for acoustic and elastic wave calculations.  Communications on Pure and Applied Mathematics, 32:313–357, 1979.
  * [12] S. Gottlieb and C. W. Shu.  Total variation diminishing Runge-Kutta schemes.  Mathematics of Computation, 67:73–85, 1998.
  * [13] G. W. Hedstrom.  Nonreflecting boundary conditions for nonlinear hyperbolic systems.  Journal of Computational Physics, 30(2):222–237, 1979.
  * [14] C. Hirsch.  Numerical Computation of Internal and External Flows.  John Wiley, 1988.
  * [15] A. Jeffrey and A. Tanuiti.  Non-linear Wave Propagation with Application to Physics and Magnetohydrodynamics.  Academic Press, 1964.
  * [16] C. Jiang, X. Feng, J. Zhang, and D. Zhong.  AMR simulations of magnetohydrodynamic problems by the CESE method in curvilinear coordinates.  Solar Physics, 267:463–491, 2010.
  * [17] G. S. Jiang and C. W Shu.  Efficient implementation of weighted ENO schemes.  Journal of Computational Physics, 126(1):202–228, 1996.
  * [18] G. S. Jiang and C. C. Wu.  A high-order WENO finite difference scheme for the equations of ideal magnetohydrodynamics.  Journal of Computational Physics, 150(2):561–594, 1999.
  * [19] L. Jiang, H. Shan, C. Liu, and M. R. Visbal.  Non-reflecting boundary conditions for DNS in curvilinear coordinates.  Recent Advances in DNS and LES. Fluid Mechanics and its Applications, 54:219–233, 1999.
  * [20] H. O. Kreiss.  Initial boundary value problems for hyperbolic systems.  Communications on Pure and Applied Mathematics, 23:277–298, 1970.
  * [21] S. Lele.  Compact finite difference schemes with spectral-like resolution.  Journal of Computational Physics, 103(1):16–42, 1992.
  * [22] B. E. Mitchell.  Direct Computation of the Sound Generated by Subsonic and Supersonic Axisymmetric Jets.  Ph.D. Dissertation, Stanford University, 1995.
  * [23] T. Poinsot and D. Veynante.  Theoretical and Numerical Combustion.  Edwards, 2005.
  * [24] T. J. Poinsot and S. K. Lele.  Boundary conditions for direct simulations of compressible viscous flows.  Journal of Computational Physics, 101:104–129, 1992.
  * [25] K. G. Powell.  An Approximate Riemann Solver for Magnetohydrodynamics. in: Hussaini M.Y., van Leer B., Van Rosendale J. (eds) Upwind and High-Resolution Schemes.  Springer, Berlin, Heidelberg, 1997.
  * [26] P. L. Roe and D. S. Balsara.  Notes on the eigensystem of magnetohydrodynamics.  SIAM Journal of Applied Mathematics, 56(1):57–67, 1996.
  * [27] C. W. Shu.  Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws.  Advanced Numerical Approximation of Nonlinear Hyperbolic Equations. Lecture Notes in Mathematics, 1697:325–423, 1998.
  * [28] G. A. Sod.  A survey of several finite difference methods for systems of nonlinear hyperbolic conservation laws.  Journal of Computational Physics, 27(1):1–31, 1978.
  * [29] K. W. Thompson.  Time dependent boundary conditions for hyperbolic systems.  Journal of Computational Physics, 68(1):1–24, 1987.
  * [30] G. Tóth.  The $\nabla\cdot{B}=0$ constraint in shock-capturing magnetohydrodynamics code.  Journal of Computational Physics, 161(2):605–652, 2000.
  * [31] A. L. Zachary and P. Colella.  A higher-order Godunov method for the equations of ideal magnetohydrodynamics.  Journal of Computational Physics, 99(2):341–347, 1992.