\section{Numerical examples}
\label{numer}
In this section, we perform numerical simulations with our positivity-preserving schemes in 1D, 2D and 3D. SSP-RK3 scheme serves as the time integrator in all the examples whereas fifth-order finite difference WENO-HCL scheme is used for solving the base MHD equations in different examples.
In multi-D, a fourth-order CT method is used to obtain a divergence-free magnetic field.
Unless otherwise stated, the gas constant is $\gamma = 5/3$ and the CFL number is $0.5$.

\subsection{Test cases in 1D}
In this subsection, we test our positivity-preserving scheme by several 1D MHD examples. 
We note that, for all the cases presented in this subsection, negative pressure or density is observed if the base MHD scheme is applied without a positivity-preserving limiter.
Here, the base MHD scheme is fifth-order WENO-HCL scheme.

\subsubsection{Vacuum shock tube test}
We first consider a 1D vacuum shock tube problem.
This example is used to demonstrate our positivity-preserving MHD solver can handle very low density and pressure. 
The initial condition is:
\begin{align}
(\rho, u_x,u_y,u_z,p,B_x,B_y,B_z) =
\begin{cases}(10^{-12},0, 0, 0, 10^{-12}, 0,0,0) & \mbox{if} \quad x<0, \\ 
( 1, 0, 0, 0, 0.5, 0,1,0) & \mbox{if} \quad x>0. \end{cases}
\end{align}
It is similar to the vacuum shock tube problem in \cite{article:Wa09}.
The computational domain is $ [-0.5,0.5]$ and zero-order extrapolation boundary conditions are used.
Shown in Figure \ref{vacuum} are the density and pressure of the solution on a mesh with $N=200$ and the highly resolved solution with $N=2000$. 
We can observe the solution of low resolution and high resolution are in good agreements. 


\subsubsection{Torsional Alfv\'en wave pulse}
We also consider the torsional Alfv\'en wave pulse problem \cite{article:BaSp99b, article:Cheng13}.
The initial condition is 
\begin{align}
(\rho, u_x,u_y,u_z,p,B_x,B_y,B_z) =
(1,10, 10 \cos \phi, 10 \sin \phi, 0.01, -10 \cos \phi, -10 \sin \phi,0),
\end{align}
where $\phi = \frac{\pi}{8} (\tanh(\frac{0.25+x}{\delta}+1))(\tanh(\frac{0.25-x}{\delta}+1))$
and $\delta = 0.005$. The computational domain is $[-0.5, 0.5]$ and periodic boundary conditions are used.
In this test problem, the initial pressure is so small 
that the problem is very sensitive to the dissipation introduced by numerical schemes. Further, the existence of a strong torsional Alfv\'en wave discontinuity makes the problem difficult to simulate. In the simulation without the  limiter, the base WENO-HCL introduced a negative pressure in a few time steps and the solutions become unphysical immediately. 
With the limiter, our scheme can simulate the problem stably and the numerical results at $t = 0.156$ are shown in Figure \ref{alfven1} and \ref{alfven2} with $N = 800$. 
Shown in the figures are plots of the energy, the thermal pressure,  $u_y$, $u_z$, $B_y$ and $B_z$. 
It is observed that our method successfully captures the two discontinuities and the results are comparable with those in \cite{article:BaSp99b, article:LiXuYa11}.
However, small bumps can still be observed around one of the discontinuities of both $u_y$ and $u_z$.
The authors in \cite{article:BaSp99b} pointed out this is because the MHD solver introduced too much numerical dissipation to keep the pressure positive.
The primary reason is the Riemann solver around the discontinuities is not selective enough.


 

\subsection{Test cases in multi-D}
In this subsection, we consider several 2D and 3D examples to demonstrate the accuracy and efficiency of our positivity-preserving multi-D MHD solver in CT framework.
In the following tests, we implement fourth-order WENO-CT2D and WENO-CT3D schemes 
as the MHD solver, to which we apply our positivity-preserving limiter.
Unless otherwise stated, we use {\bf Option 2} for the multi-D simulation. 


\subsubsection{Smooth vortex test in MHD}
We consider the smooth vortex problem with non-zero magnetic field
to demonstrate the scheme can maintain the designed accuracy within the CT framework.
We consider a modification of the smooth vortex problem considered in \cite{article:Ba04, article:LiXu12, xu2011divergence}.
The initial condition is a mean flow
\begin{align}
(\rho, u_x, u_y, u_z, p, B_x, B_y, B_z) = (1, 1, 1, 0, 1, 0, 0, 0),
\end{align}
with perturbations on $u_x$, $u_y$, $B_x$, $B_y$ and $p$:
\begin{gather*}
(\delta u_x, \delta u_y) = \frac{\kappa}{2\pi}e^{0.5(1-r^2)}(-{y},{x}), \qquad
(\delta B_x, \delta B_y) = \frac{\mu}{2\pi}e^{0.5(1-r^2)}(-{y},{x}), \\
\delta p = \frac{\mu^2(1-r^2) - \kappa^2}{8\pi^2}e^{1-r^2}.
\end{gather*}
The magnetic potential is initialized as
\begin{gather*}
A_z = \frac{\mu}{2\pi} e^{0.5(1-r^2)}
\end{gather*}
Here $r^2={x}^2 + {y}^2$. 

We set the vortex strength $\mu=5.389489439$ and $\kappa = \sqrt{2} \mu$ such that the lowest 
pressure in the center of the vortex is $5.3 \times 10^{-12}$.
Similar to \cite{article:LiXu12}, we use computational domain $(x,y) \in [-10,10]\times[-10,10]$ 
such that the error from the boundary conditions will not influence the overall convergence study.
The periodic boundary condition are used on all sides.
Because fourth-order CT steps are used, the overall scheme is fourth-order accuracy.

The $L_1$-errors and $L_\infty$-errors of the velocity and magnetic field for $t = 0.05$ are shown in Tables~\ref{tab:2dmhd}, in which one can conclude the proposed positivity-preserving scheme can maintain fourth-order accuracy as expected. 
We remark that negative pressure is observed on meshes coarser than $320 \times 320$
when the limiter is not applied.


\begin{table}
\begin{center}
\begin{Large}
    \caption{Accuracy test of the 2D vortex evolution in MHD. Shown are the $L_{1}$-errors and $L_{\infty}$-errors at time $t = 0.05$ of the density as computed by the positivity-preserving WENO-CT2D scheme at various grid resolutions. 
    The solutions converge at fourth-order accuracy.   \label{tab:2dmhd}}
    \end{Large}


\begin{tabular}{|c || c c c c|| c c c c| }
  \hline
{\normalsize{ Mesh}} & \multicolumn{4}{c||}{$u_x$} & \multicolumn{4}{c|}{$u_y$} \\ \cline{2-5} \cline{6-9}
 & {\normalsize  $L_{1}$-Error}  &  {\normalsize  Order} & {\normalsize $L_{\infty}$-Error } & {\normalsize Order }
 & {\normalsize  $L_{1}$-Error}  &  {\normalsize  Order} & {\normalsize $L_{\infty}$-Error } & {\normalsize Order }
 \\ \hline \hline
{\normalsize $40 \times 40$}    & 7.38E-04	&	-	&	1.79E-02	&	-	
								& 8.03E-04	&	-	&	1.94E-02	&	-\\
								
{\normalsize $80 \times 80$}  	&7.20E-05	&	3.35	&	4.33E-03& 	2.05	
								&7.36E-05	&	3.45	&	5.22E-03& 	1.90\\
								
{\normalsize $160 \times 160$}  &3.46E-06	&	4.38	&	1.92E-04&	4.49	
								 &3.72E-06	&	4.31	&	2.18E-04&	4.58\\
								
{\normalsize $320 \times 320$} 	&1.80E-07	&	4.27	&	1.42E-05&	3.76	
								&1.96E-07	&	4.25	&	1.64E-05&	3.73\\  \hline
 \end{tabular}
 
 \vspace{.2in}
 
\begin{tabular}{|c || c c c c|| c c c c |}
  \hline
{\normalsize{ Mesh}} & \multicolumn{4}{c||}{$B_x$} & \multicolumn{4}{c|}{$B_y$} \\ \cline{2-5} \cline{6-9}
 & {\normalsize  $L_{1}$-Error}  &  {\normalsize  Order} & {\normalsize $L_{\infty}$-Error } & {\normalsize Order }
 & {\normalsize  $L_{1}$-Error}  &  {\normalsize  Order} & {\normalsize $L_{\infty}$-Error } & {\normalsize Order }
 \\ \hline \hline
{\normalsize $40 \times 40$}    & 1.02E-03	&	-	&	1.49E-02	&	-	
								& 1.04E-03	&	-	&	1.57E-02	&	-\\
								
{\normalsize $80 \times 80$}  	&7.73E-05	&	3.73	&	1.27E-03& 	3.56	
								&7.73E-05	&	3.75	&	1.16E-03& 	3.77\\
								
{\normalsize $160 \times 160$}  &4.75E-06	&	4.03	&	8.25E-05&	3.94	
								 &4.74E-06	&	4.03	&	7.16E-05&	4.01\\
								
{\normalsize $320 \times 320$} 	&2.85E-07	&	4.06	&	7.66E-06&	3.43	
								&2.84E-07	&	4.06	&	6.36E-06&	3.49\\  \hline
 \end{tabular}
\end{center}
\end{table}

\subsubsection{Rotated vacuum shock tube problem}
\label{section:2drs}
In this example, we consider the vacuum shock tube problem rotated by an angle of $\alpha$ in a 2D domain. The initial conditions in this case are
\begin{align}
(\rho, u_{\perp},u_{\parallel},u_z,p,B_{\perp},B_{\parallel},B_z) =
\begin{cases}(10^{-12},0, 0, 0, 10^{-12}, 0,0,0) & \mbox{if} \quad \xi<0, \\ 
( 1, 0, 0, 0, 0.5, 0,1,0) & \mbox{if} \quad \xi>0. \end{cases}
\end{align}
where 
\begin{equation}
\xi = x \cos \alpha + y \sin \alpha \quad \text{and} \quad \eta = -x \sin \alpha + y \cos \alpha,
\end{equation}
where $u_{\perp}$ and $B_{\perp}$ are perpendicular to the shock interface, and $u_{\parallel}$ and $B_{\parallel}$  are parallel to the shock interface. The magnetic potential is initialized as
\begin{align}
A_z(0,\xi) =
\begin{cases} 0 & \mbox{if} \quad \xi \le 0, \\ 
- \xi & \mbox{if}  \quad \xi \ge 0. \end{cases}
\end{align}
We solve this problem by the positivity-preserving WENO-CT2D scheme on the computational domain $(x,y) \in [-0.6,0.6] \times [-0.25,0.25]$  with a $240 \times 100$ mesh. $\alpha = \tan^{-1}(0.5)$.
Zero-order extrapolation boundary conditions are used on the left and right boundaries.
On the top and bottom boundaries, all the quantities are set to describe the exact motion of the shock. 


The solutions are plotted in Figure \ref{vacuum2d}, where 1D cut of density and pressure at $y = 0$ is also plotted
to compare with the 1D highly resolved results. 
We clearly observe that the 2D solution is consistent with the 1D solution.
Without the limiter, negative density and pressure are observed in numerical solutions, 
which quickly leads to blow-up of the numerical simulation.


\subsubsection{2D blast problem}
In the blast wave problem, a strong fast magnetosonic shock formulates and propagates into the low-$\beta$  plasma background, which will likely lead to negative density or pressure in numerical solutions.
In this subsection, we first investigate a 2D version of the problem \cite{article:Bal09-ader, article:BaSp99a, article:LiXuYa11}.
The computational domain is $(x,y) \in [-0.5,0.5] \times [-0.5, 0.5]$ with outflow boundary conditions on all the four sides. 

The initial conditions of the problem consist of an initial background:
\begin{equation}
(\rho, u_x, u_y, u_z, p, B_x, B_y, B_z) = (1, 0, 0, 0, 0.1, 100/\sqrt{2 \pi},100/\sqrt{2 \pi},0),
\end{equation}
and a circular pressure pulse $p = 1000$ within a radius $r=0.1$ from the center of the domain. The initial scalar magnetic potential is simply given by
\begin{align}
A_z =  100/\sqrt{2 \pi} y - 100/\sqrt{2 \pi} x.
\end{align}
The solution is computed on a $256 \times 256$ mesh. Shown in Figure~\ref{blast2} are plots of the solutions. 
The solution shows good agreement with those in \cite{article:Bal09-ader,article:LiXuYa11}.

In Table~\ref{tab:comp}, we use this example to compare four different schemes, WENO-HCL, WENO-CT-OP1, WENO-CT-OP2 and PP-WENO-CT-OP2. 
Here WENO-HCL is referred to the base WENO-HCL scheme without CT or the limiter.
WENO-CT-OP1 is referred to the WENO-CT2D scheme choosing {\bf Option 1} without the limiter. 
WENO-CT-OP2 is referred to the WENO-CT2D scheme choosing {\bf Option 2} without the limiter. 
Finally, PP-WENO-CT-OP2 is referred to the positivity-preserving WENO-CT2D scheme with {\bf Option 2} chosen.

From Table~\ref{tab:comp}, we observe that the base WENO-HCL scheme is unstable for the resolution $150 \times 150$ and becomes stable in the higher resolutions.
WENO-CT-OP1 is unstable for each resolution and applying the positivity-preserving limiter will not be able to stabilize this because the negative pressure is from the correction step of the magnetic field.
WENO-CT-OP2 is stable in the lower resolution but become unstable for the resolution $256 \times 256$, and negative pressure is observed in all the resolutions.
Finally, the positivity-preserving WENO-CT scheme is stable for all the resolutions.
From those results, it is very clear that the positivity-preserving WENO-CT scheme is the most stable methods in those four methods.

Another concern in the CT framework is that the energy is not conservative in our positivity-preserving WENO-CT scheme due to {\bf Option 2}. 
We also use this example to study this issue. 
We compare the results by the base WENO-HCL scheme and the 
positivity-preserving WENO-HCL scheme.
In Figure~\ref{blast_hcl}, we show the results by the base WENO-HCL scheme with the same resolution as Figure~\ref{blast2}.
The results look similar to those by the positivity-preserving WENO-CT scheme,
except there are some unphysical oscillations around the center region in Figure~\ref{blast_hcl}.
That is due to the divergence error in the base scheme. 
If we plot the divergence error in the time domain, 
we can clearly see the divergence error of the positivity-preserving WENO-CT scheme
stays around $10^{-12}$, while the error of the WENO-HCL scheme is around $10^0$ during the simulation.
Here the divergence error is defined as $L_1$-norm of $\divg \B$,
where the numerical $\divg$ operator is defined as a regular fourth-order central finite difference discretization.
As a common drawback in CT framework when {\bf Option 2} is chosen,
the correction step leads to a loss of the conservation of the total energy.
The results are plotted in Figure~\ref{blast_hcl}, where we can see 
the relative total energy error is around $10^{-3}$ 
while the conservative WENO-HCL scheme has an error about $10^{-12}$.
But we find this loss of total energy will decrease as the mesh is refined.
However, we remark that the conservation of the total energy is important for some problems, 
such as those involving nonlinear strong discontinuities.
A high-order positivity-preserving conservative scheme with  
the divergence error controlled will be part of our future work.
However, it is very difficulty, if not impossible, in the CT framework to satisfy 
all the requirements simultaneously.
A better way to control the divergence error is needed for this purpose.

\begin{table}
\begin{center}
\begin{Large}
\caption{Comparisons of different schemes solving the 2D blast problem. 
The column of ``Positivity'' lists if a negative solution is observed in the simulation. The column of ``Stability'' lists if the simulation run stably to $t = 0.01$.
In order to make a fair comparison, the positivity of density and pressure is only checked at each time step $t^n$. 
\label{tab:comp}}
\end{Large}
\begin{tabular}{|c || c c || c c|| c c || c c |}
  \hline
{\normalsize{ Mesh}} & \multicolumn{2}{c||}{WENO-HCL} & \multicolumn{2}{c||}{WENO-CT-OP1} & \multicolumn{2}{c||}{WENO-CT-OP2} & \multicolumn{2}{c|}{PP-WENO-CT-OP2} \\ 
\cline{2-5} \cline{6-9}
 & {\normalsize  Positivity}  &  {\normalsize  Stability} & {\normalsize Positivity } & {\normalsize Stability }
 & {\normalsize  Positivity}  &  {\normalsize  Stability} & {\normalsize Positivity } & {\normalsize Stability }
 \\ \hline \hline
{\normalsize $150 \times 150$}  & No	&	No	&	No	&	No
								& No	&	Yes	&	Yes	&	Yes\\
								
{\normalsize $200 \times 200$}  & No	&	Yes	&	No	& 	No
								& No	&	Yes	&	Yes & 	Yes\\
								
{\normalsize $256 \times 256$}  &No	&	Yes	&	No	&	No	
								&No	&	No	&	Yes &	Yes\\
 \hline
 \end{tabular}
 \end{center}
 \end{table}


\subsubsection{3D blast problem}
The last problem we investigate is a fully 3D version of the blast problem.
It is used to test the behavior of the positivity-preserving WENO-CT3D scheme.
The initial conditions consist of an initial background:
\begin{equation}
(\rho, u_x, u_y, u_z, p, B_x, B_y, B_z) = (1, 0, 0, 0, 0.1, 100/\sqrt{4 \pi}/\sqrt{2},100/\sqrt{4 \pi}/\sqrt{2},0)
\end{equation}
and a spherical pressure pulse $p = 1000$ within a radius $r=0.1$ from the centre of the domain.
The initial conditions for the magnetic potential are
\begin{align}
\Av(0,x,y,z) = (0, 0, 100/\sqrt{4 \pi}/\sqrt{2} y - 100/\sqrt{4 \pi}/\sqrt{2} x ).
\end{align}
The computational domain is $[-0.5,0.5]^3$. 
Outflow boundary conditions are used on all the sides. The numerical simulation is performed on a $150 \times 150 \times 150$ mesh. 
In Figure \ref{blast3d} we show the results of the solutions cut at $z = 0$.
To distinguish this 3D case from the 2D blast case, 
we also present the 3D plots of the density and pressure in Figure \ref{3dplot}, 
which clearly indicates its spherical structures.
The solution is comparable to the 3D results in \cite{article:GaSt08, article:MiTz10, article:Zi04}.
We note that negative pressure is observed at time $t = 0.0033$ if the positivity preserving limiter is not applied.


\begin{figure}
\begin{center}
\begin{tabular}{cc}
  	(a)\includegraphics[width=0.42\textwidth]{dennew-eps-converted-to} & 
	(b)\includegraphics[width=0.42\textwidth]{prenew-eps-converted-to} \\
\end{tabular}
      \caption{Vacuum shock tube problem. 
      Shown in these panels are plots at time $t=0.1$ of 
      (a) the density and
      (b) the thermal pressure. 
      The blue circle is a solution solved on a mesh with $N = 200$.
      The solid line is a highly resolved solution with $N = 2000$.     
      \label{vacuum}}
\end{center}
\end{figure}

\begin{figure}
\begin{center}
\begin{tabular}{cc}
  	(a)\includegraphics[width=0.42\textwidth]{eng-eps-converted-to} & 
	(b)\includegraphics[width=0.42\textwidth]{pre-eps-converted-to}
\end{tabular}
      \caption{The torsional Alfv\'en wave pulse. 
      Shown in these panels are plots at time $t=0.156$ of 
      (a) the energy and
      (b) the thermal pressure.
      The solution was obtained on a mesh with $N=800$.     
      \label{alfven1}}
\end{center}
\end{figure}

\begin{figure}
\begin{center}
\begin{tabular}{cc}
  	(a)\includegraphics[width=0.42\textwidth]{u2-eps-converted-to} & 
	(b)\includegraphics[width=0.42\textwidth]{u3-eps-converted-to} \\
	(c)\includegraphics[width=0.42\textwidth]{b2-eps-converted-to} & 
	(d)\includegraphics[width=0.42\textwidth]{b3-eps-converted-to}
\end{tabular}
      \caption{The torsional Alfv\'en wave pulse. 
      Shown in these panels are plots at time $t=0.156$ of 
      (a) $u_y$,
      (b) $u_z$,
      (c) $B_y$ and
      (d) $B_z$. 
      The solution was obtained on a mesh with $N=800$.    
      \label{alfven2}}
\end{center}
\end{figure}

\begin{figure}
\begin{center}
\begin{tabular}{cc}
  	(a)\includegraphics[width=0.42\textwidth]{den2d-eps-converted-to} & 
	(b)\includegraphics[width=0.42\textwidth]{dencut-eps-converted-to} \\
	(c)\includegraphics[width=0.42\textwidth]{pre2d-eps-converted-to} & 
	(d)\includegraphics[width=0.42\textwidth]{precut-eps-converted-to}
\end{tabular}
      \caption{Rotated vacuum shock tube problem. 
      Shown in these panels are plots at time $t=0.1$ of 
      (a) the density,
      (b) the density cut at y = 0,
      (c) the pressure and
      (d) the pressure cut at y = 0. 
      40 equally spaced contours are used for the contour plots. The solid lines in (b) and (d)
      are 1D highly resolved solutions.
      The solution was obtained on a $240 \times 100$ mesh.     
      \label{vacuum2d}}
\end{center}
\end{figure}


\begin{figure}
\begin{center}
\begin{tabular}{cc}
  	(a)\includegraphics[width=0.42\textwidth]{blast2den} & 
	(b)\includegraphics[width=0.42\textwidth]{blast2pre} \\
	(c)\includegraphics[width=0.42\textwidth]{blast2unm} & 
	(d)\includegraphics[width=0.42\textwidth]{blast2bnm}
\end{tabular}
      \caption{2D blast problem. 
      Shown in these panels are plots at time $t=0.01$ of 
      (a) the density,
      (b) the thermal pressure,
      (c) the norm of velocity and
      (d) the magnetic pressure. 
      40 equally spaced contours are used for each plot.
      The solution was obtained on a $256 \times 256$ mesh by positivity-preserving WENO-CT scheme.    
      \label{blast2}}
\end{center}
\end{figure}


\begin{figure}
\begin{center}
\begin{tabular}{cc}
  	(a)\includegraphics[width=0.42\textwidth]{blast_hclden} & 
	(b)\includegraphics[width=0.42\textwidth]{blast_hclpre} \\
	(c)\includegraphics[width=0.42\textwidth]{blast_hclunm} & 
	(d)\includegraphics[width=0.42\textwidth]{blast_hclbnm}
\end{tabular}
      \caption{2D blast problem. 
      Shown in these panels are plots at time $t=0.01$ of 
      (a) the density,
      (b) the thermal pressure,
      (c) the norm of velocity and
      (d) the magnetic pressure. 
      40 equally spaced contours are used for each plot.
      The solution was obtained on a $256 \times 256$ mesh by WENO-HCL scheme.    
      \label{blast_hcl}}
\end{center}
\end{figure}

\begin{figure}
\begin{center}
\begin{tabular}{cc}
  	(a)\includegraphics[width=0.42\textwidth]{divgerr-eps-converted-to} & 
	(b)\includegraphics[width=0.42\textwidth]{engerr-eps-converted-to}
\end{tabular}
      \caption{Comparisons between WENO-HCL scheme and positivity-preserving WENO-CT scheme for the 2D blast problem. 
      Shown in these panels are plots of 
      (a) the divergence error and 
      (b) the relative total energy error in the time domain $t \in [0,0.01]$.     
      \label{blast_comp}}
\end{center}
\end{figure}


\begin{figure}
\begin{center}
\begin{tabular}{cc}
  	(a)\includegraphics[width=0.42\textwidth]{blast3den} & 
	(b)\includegraphics[width=0.42\textwidth]{blast3pre} \\
	(c)\includegraphics[width=0.42\textwidth]{blast3unm} & 
	(d)\includegraphics[width=0.42\textwidth]{blast3bnm}
\end{tabular}
      \caption{3D blast problem. 
      Shown in these panels are plots at time $t=0.01$ and cut at $z = 0$ of 
      (a) the density,
      (b) the thermal pressure,
      (c) the norm of velocity and
      (d) the magnetic pressure. 
      40 equally spaced contours are used for each plot.
      The solution was obtained on a $150 \times 150 \times 150$ mesh.     
      \label{blast3d}}
\end{center}
\end{figure}

\begin{figure}
\begin{center}
\begin{tabular}{cc}
(a) \includegraphics[width=0.42\textwidth]{den3d} &
(b) \includegraphics[width=0.42\textwidth]{pre3d}
\end{tabular}
\caption{3D blast problem. 
      Shown in these panels are 3D plots at time $t=0.01$ of 
      (a) the density and
      (b) the thermal pressure.
      The solution was obtained on a $150 \times 150 \times 150$ mesh.     
      \label{3dplot}
}
\end{center}
\end{figure}
{"config": "arxiv", "file": "1406.5098/numerical.tex"}