A New Locally Divergence-Free Path-Conservative Central-Upwind Scheme for Ideal and Shallow Water Magnetohydrodynamics

Alina Chertock, Alexander Kurganov, Michael Redle, Kailiang Wu Department of Mathematics and Center for Research in Scientific Computing, North Carolina State University, Raleigh, NC 27695, USA; [email protected]Department of Mathematics, SUSTech International Center for Mathematics, and Guangdong Provincial Key Laboratory of Computational Science and Material Design, Southern University of Science and Technology, Shenzhen 518055, China; [email protected]Department of Mathematics, North Carolina State University, Raleigh, NC 27695, USA; [email protected]Department of Mathematics and SUSTech International Center for Mathematics, Southern University of Science and Technology, and National Center for Applied Mathematics Shenzhen (NCAMS), Shenzhen 518055, China; [email protected]

Abstract

We develop a new second-order unstaggered path-conservative central-upwind (PCCU) scheme for ideal and shallow water magnetohydrodynamics (MHD) equations. The new scheme possesses several important properties: it locally preserves the divergence-free constraint, it does not rely on any (approximate) Riemann problem solver, and it robustly produces high-resolution and non-oscillatory results. The derivation of the scheme is based on the Godunov-Powell nonconservative modifications of the studied MHD systems. The local divergence-free property is enforced by augmenting the modified systems with the evolution equations for the corresponding derivatives of the magnetic field components. These derivatives are then used to design a special piecewise linear reconstruction of the magnetic field, which guarantees a non-oscillatory nature of the resulting scheme. In addition, the proposed PCCU discretization accounts for the jump of the nonconservative product terms across cell interfaces, thereby ensuring stability. We test the proposed PCCU scheme on several benchmarks for both ideal and shallow water MHD systems. The obtained numerical results illustrate the performance of the new scheme, its robustness, and its ability not only to achieve high resolution, but also preserve the positivity of computed quantities such as density, pressure, and water depth.

Keywords: Ideal magnetohydrodynamics, shallow water magnetohydrodynamics, divergence-free constraints, path-conservative central-upwind scheme, nonconservative hyperbolic systems of nonlinear PDEs.

AMS subject classification: 65M08, 76W05, 76M12, 86-08, 35L65.

1 Introduction

This paper focuses on developing a novel numerical method for magnetohydrodynamic (MHD) systems, widely used in many applications, such as astrophysics, plasma physics, space physics, and engineering. In these models, fluid dynamics equations are coupled with the equations for the magnetic field, which satisfies the divergence-free condition – a physically-exact constraint, that is, if initially, the divergence of the magnetic field is zero, then it must remain zero for all times. When deriving numerical methods for MHD systems, the divergence-free condition must be handled with care, as neglecting an identically-zero divergence on a discrete level may lead to severe numerical instabilities and/or nonphysical structures in the numerical solution; see, e.g., [5, 6, 38, 56]. In addition, like other hyperbolic systems of conservation and balance laws, the MHD systems typically develop very complicated nonsmooth solution structures containing shock waves, rarefactions, and contact discontinuities, as well as their interactions.

In the past decades, various numerical techniques have been developed to deal with the divergence-free constraint for MHD systems. An early effort in this direction is the projection method [6], which is a post-processing divergence correction procedure that uses Hodge decomposition to project the non-divergence-free magnetic field into a divergence-free subspace by solving an elliptic Poisson equation. Another widely used approach is the constrained transport (CT) method, which was proposed in [21] for simulating MHD flows. This method preserves a specific discrete version of divergence-free condition on staggered grids, and its variants were further developed by researchers within various frameworks; see, e.g., [17, 12, 54, 5, 44, 24, 64]. Unstaggered CT methods were also developed (see, e.g., [53, 29, 45, 30, 10]), and they are usually based on numerically evolving the magnetic potential and computing the divergence-free magnetic field through the (discrete) curl of the magnetic potential. In addition, locally divergence-free discontinuous Galerkin methods that enforce the zero divergence of the magnetic field within each cell were developed in [38, 65]. In recent years, globally divergence-free high-order methods were also proposed to enforce the exact zero divergence of the magnetic field within the finite-volume or (central) discontinuous Galerkin framework; see, e.g., [2, 4, 19, 40, 39, 22, 3].

There is also a different class of schemes that reduce the divergence errors but do not explicitly enforce any divergence-free constraint. In the context of the ideal MHD equations, these methods, typically referred to as eight-wave methods, were proposed by Powell [51, 49, 50] based on a proper discretization of a modified, nonconservative ideal MHD model. This model was first introduced by Godunov [26] for entropy symmetrization. Compared to the conservative ideal MHD equations, the modified model contains extra nonconservative source terms (referred to as Godunov-Powell source terms in the following), which are proportional to the divergence of the magnetic field. These source terms change the character of the MHD equations, making the modified model Galilean invariant [14], symmetrizable [26], and helpful in designing entropy stable schemes (see, e.g., [9, 16, 43]). In [49], it was also noticed that the conservative ideal MHD equations are weakly hyperbolic, and thus such source terms should be added to recover the missing eigenvector. As demonstrated in [51, 50], the inclusion of the source terms ensures that the magnetic divergence is advected with the flow, and the numerical divergence errors are also expected to be advected and would not accumulate. This makes the eight-wave methods capable of controlling the divergence error, although certain drawbacks may arise due to a nonconservative nature of the Godunov-Powell modified ideal MHD equations; see [56]. As recently discovered in [60], a discrete divergence-free condition is closely related to the positivity-preserving property of numerical schemes for the ideal MHD equations. Furthermore, locally divergence-free positivity-preserving schemes [62, 63] for the Godunov-Powell modified ideal MHD model can be obtained via geometric quasilinearization [61]. Another class of divergence-controlling schemes is the so-called hyperbolic divergence-cleaning method [13], which introduces a mixed hyperbolic-parabolic equation to damp the divergence errors away instead of enforcing an exactly divergence-free magnetic field.

A variant of the MHD equations—known as the shallow water MHD system—has also become a model of significant numerical interest over the last few decades; see, e.g., [37, 67, 48] and references therein. First introduced in the context of a solar tachocline in [25] and now used in several astrophysical and geophysical contexts, this variant is fully derived from the ideal MHD equations under the assumptions of constant density and magneto-hydrostatic equilibrium; see [15, 66]. An assortment of numerical methods have additionally been explored to treat divergence errors for this system—such as space-time conservation element solution element methods (CE/SE) in [1, 52], an evolution Galerkin scheme in [32], and entropy-stable schemes in [59, 18], to name a few.

Despite these advances, devising highly accurate, stable, and at the same time, robust numerical methods capable of preserving the divergence-free condition at a discrete level is still a challenging task. Our main goal is to develop such a scheme. To this end, we consider the Godunov-Powell modified ideal and shallow water MHD models and supplemented them with additional equations obtained by differentiating the magnetic field equations in space: the latter will help to ensure local divergence-free conditions. The resulting augmented MHD systems will be nonconservative and rather complicated to be solved by an upwind numerical method, that is, by a method relying on a solution of (generalized) Riemann problems. Instead, we numerically solve the augmented MHD systems using second-order unstaggered finite-volume Riemann-problem-solver-free path-conservative central-upwind (PCCU) schemes, which were introduced in [8] as a black-box solver for nonconservative hyperbolic systems of PDEs. PCCU schemes are a path-conservative extension of the central-upwind (CU) schemes, which were developed in [35, 34, 33] for general multidimensional hyperbolic systems of conservation laws. We enforce the local divergence-free condition with the help of a special piecewise linear reconstruction of the magnetic field variables. The resulting scheme produces highly accurate and non-oscillatory results for ideal and shallow water MHD systems.

The paper consists of two parts: In §2, we study the ideal MHD equations, while §3 is devoted to the shallow water MHD system. The Godunov-Powell modifications and the augmented forms of the studied systems are presented in §2.1 and §3.1, the new numerical methods for the resulting augmented systems are introduced in §2.2 and §3.2, and the results of the preformed numerical experiments are reported in §2.3 and §3.3. We conclude the paper and discuss some of our future research plans in §4.

2 Ideal MHD Equations

2.1 Governing Equations

The ideal MHD equations read as

	$\displaystyle\rho_{t}+\nabla\cdot{\left(\rho\bm{u}\right)}=0,$
	$\displaystyle(\rho\bm{u})_{t}+\nabla\cdot\Big{[}\rho\bm{u}\bm{u}^{\top}+\Big{(}p+\frac{1}{2}{\left\|\bm{b}\right\|}^{2}\Big{)}I-\bm{b}\bm{b}^{\top}\Big{]}=\bm{0},$
	$\displaystyle\bm{b}_{t}-\nabla\times(\bm{u}\times\bm{b})=\bm{0},$
	$\displaystyle{\cal E}_{t}+\nabla\cdot\Big{[}\Big{(}{\cal E}+p+\frac{1}{2}{\left\|\bm{b}\right\|}^{2}\Big{)}\bm{u}-\bm{b}(\bm{u}\cdot\bm{b})\Big{]}=0,$

where $t$ represents the time, $\rho$ is the density, $p$ is the pressure, $\bm{u}=(u,v,w)^{\top}$ represents the fluid velocity, $\bm{b}=(b_{1},b_{2},b_{3})^{\top}$ is the magnetic field, and ${\cal E}$ is the total energy. Additionally, $I$ is the $3\times 3$ identity matrix, $\gamma$ represents the ratio of specific heats, and finally, the system is completed through the equation of state (EOS)

{\cal E}=\frac{p}{\gamma-1}+\frac{\rho}{2}{\left|\bm{u}\right|}^{2}+\frac{1}{2}{\left|\bm{b}\right|}^{2}.

(2.2)

where $|\cdot|$ represents the Euclidean norm. It is easy to show that

\nabla\cdot\bm{b}=0,

(2.3)

as long as initially the magnetic field is divergence-free.

As mentioned in §1, we follow a commonly-used approach and, instead of considering (LABEL:2.1f), we develop a new numerical method for the Godunov-Powell modified ideal MHD equations:

		$\displaystyle\rho_{t}+\nabla\cdot{\left(\rho\bm{u}\right)}=0,$		(2.4)
		$\displaystyle(\rho\bm{u})_{t}+\nabla\cdot\Big{[}\rho\bm{u}\bm{u}^{\top}+\Big{(}p+\frac{1}{2}{\left\|\bm{b}\right\|}^{2}\Big{)}I-\bm{b}\bm{b}^{\top}\Big{]}=-\bm{b}{\left(\nabla\cdot\bm{b}\right)},$
		$\displaystyle\bm{b}_{t}-\nabla\times(\bm{u}\times\bm{b})=-\bm{u}{\left(\nabla\cdot\bm{b}\right)},$
		$\displaystyle{\cal E}_{t}+\nabla\cdot\Big{[}\Big{(}{\cal E}+p+\frac{1}{2}{\left\|\bm{b}\right\|}^{2}\Big{)}\bm{u}-\bm{b}(\bm{u}\cdot\bm{b})\Big{]}=-{\left(\bm{u}\cdot\bm{b}\right)}{\left(\nabla\cdot\bm{b}\right)},$

which is closed with the help of the same EOS (2.2). Note that, theoretically, the Godunov-Powell source terms $-\bm{b}{\left(\nabla\cdot\bm{b}\right)}$ , $-\bm{u}{\left(\nabla\cdot\bm{b}\right)}$ , and $-{\left(\bm{u}\cdot\bm{b}\right)}{\left(\nabla\cdot\bm{b}\right)}$ on the right hand side (RHS) of (2.4), are zero due to the divergence-free condition (2.3). However, when numerically solving (2.4) with EOS (2.2), these added relaxation terms help to reduce the divergence errors and enhance the robustness; see, e.g., [51, 50, 23, 31, 58, 62, 63]. It is worth noting that, although our proposed schemes are locally divergence-free, there are jumps of normal magnetic component across cell interfaces, and the inclusion of these extra source terms can help to control the (weak) divergence errors at cell interfaces.

In this paper, we restrict our attention to the 2-D case, where all the quantities of interest depend on the spatial variables $x$ and $y$ and time $t$ only. In this case, the divergence-free condition (2.3) reads as $(b_{1})_{x}+(b_{2})_{y}=0$ and one of the goals in the development of a good numerical method for the ideal MHD system (2.4), (2.2) is to enforce this condition at the discrete level. In order to achieve this goal, we introduce the new variables $A:=(b_{1})_{x}$ and $B:=(b_{2})_{y}$ , and differentiate the induction equation for $b_{1}$ and $b_{2}$ in (2.4) with respect to $x$ and $y$ , respectively, to obtain the following two evolution equations for $A$ and $B$ :

		$\displaystyle A_{t}+\big{(}uA-b_{2}u_{y}\big{)}_{x}+\big{(}vA+b_{1}v_{x}\big{)}_{y}=0,$		(2.5)
		$\displaystyle B_{t}+\big{(}uB+b_{2}u_{y}\big{)}_{x}+\big{(}vB-b_{1}v_{x}\big{)}_{y}=0.$		(2.5)

From now on, we will add these equations to the Godunov-Powell modified ideal MHD equations and will numerically solve the augmented system (2.4)–(2.5), (2.2). Even though the number of equations to be discretized has been increased, adding the equations in (2.5) makes it easier to control the divergence-free constraint, which now reads as $A+B=0$ .

Before introducing the numerical method, we write the augmented ideal MHD system (2.4)–(2.5) in the following vector form:

\bm{U}_{t}+\bm{F}(\bm{U})_{x}+\bm{G}(\bm{U})_{y}=Q(\bm{U})\bm{U}_{x}+R(\bm{U})\bm{U}_{y},

(2.6)

where $\bm{U}=(\rho,\rho u,\rho v,\rho w,b_{1},b_{2},b_{3},{\mathcal{E}},A,B)^{\top}$ ,

		$\displaystyle\bm{F}(\bm{U})=\begin{pmatrix}\rho u\\ \rho u^{2}+p+\frac{1}{2}{\left\|\bm{b}\right\|}^{2}-b_{1}^{2}\\ \rho uv-b_{1}b_{2}\\ \rho uw-b_{1}b_{3}\\ 0\\ ub_{2}-vb_{1}\\ ub_{3}-wb_{1}\\ \big{(}{\cal E}+p+\frac{1}{2}{\left\|\bm{b}\right\|}^{2}\big{)}u-(\bm{u}\cdot\bm{b})b_{1}\\ uA-b_{2}u_{y}\\ uB+b_{2}u_{y}\end{pmatrix},\leavevmode\nobreak\ \leavevmode\nobreak\ Q(\bm{U})=\begin{pmatrix}0&0&0&0&0&0&0&0&0&0\\ 0&0&0&0&-b_{1}&0&0&0&0&0\\ 0&0&0&0&-b_{2}&0&0&0&0&0\\ 0&0&0&0&-b_{3}&0&0&0&0&0\\ 0&0&0&0&-u&0&0&0&0&0\\ 0&0&0&0&-v&0&0&0&0&0\\ 0&0&0&0&-w&0&0&0&0&0\\ 0&0&0&0&-\bm{u}\cdot\bm{b}&0&0&0&0&0\\ 0&0&0&0&0&0&0&0&0&0\\ 0&0&0&0&0&0&0&0&0&0\end{pmatrix},$
		$\displaystyle\bm{G}(\bm{U})=\begin{pmatrix}\rho v\\ \rho uv-b_{1}b_{2}\\ \rho v^{2}+p+\frac{1}{2}{\left\|\bm{b}\right\|}^{2}-b_{2}^{2}\\ \rho vw-b_{2}b_{3}\\ vb_{1}-ub_{2}\\ 0\\ vb_{3}-wb_{2}\\ \big{(}{\cal E}+p+\frac{1}{2}{\left\|\bm{b}\right\|}^{2}\big{)}v-(\bm{u}\cdot\bm{b})b_{2}\\ vA+b_{1}v_{x}\\ vB-b_{1}v_{x},\end{pmatrix},\leavevmode\nobreak\ \leavevmode\nobreak\ R(\bm{U})=\begin{pmatrix}0&0&0&0&0&0&0&0&0&0\\ 0&0&0&0&0&-b_{1}&0&0&0&0\\ 0&0&0&0&0&-b_{2}&0&0&0&0\\ 0&0&0&0&0&-b_{3}&0&0&0&0\\ 0&0&0&0&0&-u&0&0&0&0\\ 0&0&0&0&0&-v&0&0&0&0\\ 0&0&0&0&0&-w&0&0&0&0\\ 0&0&0&0&0&-\bm{u}\cdot\bm{b}&0&0&0&0\\ 0&0&0&0&0&0&0&0&0&0\\ 0&0&0&0&0&0&0&0&0&0\end{pmatrix}.$

2.2 Numerical Method

We introduce a uniform Cartesian mesh consisting of the finite-volume cells $C_{j,k}=[x_{j-\frac{1}{2}},x_{j+\frac{1}{2}}]\times[y_{k-\frac{1}{2}},y_{k+\frac{1}{2}}]$ with $x_{j+\frac{1}{2}}-x_{j-\frac{1}{2}}\equiv\Delta x$ and $y_{k+\frac{1}{2}}-y_{k-\frac{1}{2}}\equiv\Delta y$ . We assume that at a certain time level $t$ , the computed solution realized in terms of its cell averages,

\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}_{j,k}\approx\frac{1}{\Delta x\Delta y}\iint\limits_{C_{j,k}}\bm{U}(x,y,t)\,{\rm d}x\,{\rm d}y,

is available. Notice that the dependence of $\,\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}_{j,k}$ and many other indexed quantities on $t$ is omitted here and throughout the rest of the paper for the sake of brevity.

The cell averages $\,\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}_{j,k}$ are evolved in time by implementing a dimension-by-dimension extension of the PCCU scheme from [8], which results in the following semi-discretization of (2.6)–(LABEL:2.10):

$\displaystyle\frac{{\rm d}}{{\rm d}t}\,\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}_{j,k}=$	$\displaystyle-\frac{1}{\Delta x}\bigg{[}\bm{{\cal F}}_{{j+\frac{1}{2}},k}-\bm{{\cal F}}_{{j-\frac{1}{2}},k}-\bm{Q}_{j,k}$
	$\displaystyle\hskip 36.98866pt-\frac{s_{{j-\frac{1}{2}},k}^{+}}{s_{{j-\frac{1}{2}},k}^{+}-s_{{j-\frac{1}{2}},k}^{-}}\,\bm{Q}_{\bm{\Psi},{j-\frac{1}{2}},k}+\frac{s_{{j+\frac{1}{2}},k}^{-}}{s_{{j+\frac{1}{2}},k}^{+}-s_{{j+\frac{1}{2}},k}^{-}}\,\bm{Q}_{\bm{\Psi},{j+\frac{1}{2}},k}\bigg{]}$
		(2.8)
	$\displaystyle-\frac{1}{\Delta y}\bigg{[}\bm{{\cal G}}_{j,{k+\frac{1}{2}}}-\bm{{\cal G}}_{j,{k-\frac{1}{2}}}-\bm{R}_{j,k}$
	$\displaystyle\hskip 36.98866pt-\frac{s_{j,{k-\frac{1}{2}}}^{+}}{s_{j,{k-\frac{1}{2}}}^{+}-s_{j,{k-\frac{1}{2}}}^{-}}\,\bm{R}_{\bm{\Psi},j,{k-\frac{1}{2}}}+\frac{s_{j,{k+\frac{1}{2}}}^{-}}{s_{j,{k+\frac{1}{2}}}^{+}-s_{j,{k+\frac{1}{2}}}^{-}}\,\bm{R}_{\bm{\Psi},j,{k+\frac{1}{2}}}\bigg{]}.$

Here,

	$\displaystyle\bm{{\cal F}}_{{j+\frac{1}{2}},k}$	$\displaystyle=\frac{s_{{j+\frac{1}{2}},k}^{+}\bm{F}\big{(}\bm{U}^{\rm E}_{j,k}\big{)}-s_{{j+\frac{1}{2}},k}^{-}\bm{F}\big{(}\bm{U}^{\rm W}_{j+1,k}\big{)}}{s_{{j+\frac{1}{2}},k}^{+}-s_{{j+\frac{1}{2}},k}^{-}}+\frac{s_{{j+\frac{1}{2}},k}^{+}s_{{j+\frac{1}{2}},k}^{-}}{s_{{j+\frac{1}{2}},k}^{+}-s_{{j+\frac{1}{2}},k}^{-}}\left(\bm{U}^{\rm W}_{j+1,k}-\bm{U}^{\rm E}_{j,k}\right),$		(2.9)
	$\displaystyle\bm{{\cal G}}_{j,{k+\frac{1}{2}}}$	$\displaystyle=\frac{s_{j,{k+\frac{1}{2}}}^{+}\bm{G}\big{(}\bm{U}^{\rm N}_{j,k}\big{)}-s_{j,{k+\frac{1}{2}}}^{-}\bm{G}\big{(}\bm{U}^{\rm S}_{j,k+1}\big{)}}{s_{j,{k+\frac{1}{2}}}^{+}-s_{j,{k+\frac{1}{2}}}^{-}}+\frac{s_{j,{k+\frac{1}{2}}}^{+}s_{j,{k+\frac{1}{2}}}^{-}}{s_{j,{k+\frac{1}{2}}}^{+}-s_{j,{k+\frac{1}{2}}}^{-}}\left(\bm{U}^{\rm S}_{j,k+1}-\bm{U}^{\rm N}_{j,k}\right)$		(2.9)

are the CU numerical fluxes from [36], $\bm{U}^{\rm E,\rm W,\rm N,\rm S}_{j,k}$ are the reconstructed point values at the cell interfaces of cell $C_{j,k}$ (see §2.2.1 for details), $s_{{j+\frac{1}{2}},k}^{\pm}$ and $s_{j,{k+\frac{1}{2}}}^{\pm}$ are the one-sided local speeds of propagation in the $x$ - and $y$ -directions, respectively (see §2.2.2 for details), and $\bm{Q}_{j,k}$ , $\bm{R}_{j,k}$ , $\bm{Q}_{\bm{\Psi},{j+\frac{1}{2}},k}$ , and $\bm{R}_{\bm{\Psi},j,{k+\frac{1}{2}}}$ denote the discretizations of the nonconservative products on the RHS of (2.6) (see §2.2.3 for details).

We point out that (2.8) is a system of ODEs, which should be numerically integrated in time by an appropriate ODE solver. In the numerical experiments reported in §2.3, we have used the explicit three-stage third-order strong stability preserving (SSP) Runge-Kutta method; see, e.g., [28, 27].

2.2.1 Piecewise Linear Reconstruction

Equipped with the cell averages $\,\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}_{j,k}$ , we first use the EOS (2.2) and compute the approximate point values of $u$ , $v$ , $w$ , and $p$ at the cell centers:

		$\displaystyle u_{j,k}=\frac{(\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\rho u$\kern 0.0pt}}})_{j,k}}{\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}_{j,k}},\quad v_{j,k}=\frac{(\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\rho v$\kern 0.0pt}}})_{j,k}}{\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}_{j,k}},\quad w_{j,k}=\frac{(\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\rho w$\kern 0.0pt}}})_{j,k}}{\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}_{j,k}},$
		$\displaystyle p_{j,k}={\left(\gamma-1\right)}\bigg{[}\,\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\cal E$\kern 0.0pt}}}_{j,k}-\frac{1}{2}\,\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}_{j,k}\left(u_{j,k}^{2}+v_{j,k}^{2}+w_{j,k}^{2}\right)-\frac{1}{2}\Big{(}(\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$b_{1}$\kern 0.0pt}}})_{j,k}^{2}+(\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$b_{2}$\kern 0.0pt}}})_{j,k}^{2}+(\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$b_{3}$\kern 0.0pt}}})_{j,k}^{2}\Big{)}\bigg{]}.$

We then introduce a new set of discrete variables.

\bm{W}_{j,k}:=(\,\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\rho$\kern 0.0pt}}}_{j,k},u_{j,k},v_{j,k},w_{j,k},(\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$b_{1}$\kern 0.0pt}}})_{j,k},(\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$b_{2}$\kern 0.0pt}}})_{j,k},(\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$b_{3}$\kern 0.0pt}}})_{j,k},p_{j,k},\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$A$\kern 0.0pt}}}_{j,k},\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$B$\kern 0.0pt}}}_{j,k})^{\top},

and compute the cell interface point values $\bm{W}^{\rm E,\rm W,\rm N,\rm S}_{j,k}$ using a proper conservative piecewise linear reconstruction

\widetilde{\bm{W}}(x,y)=\bm{W}_{j,k}+(\bm{W}_{x})_{j,k}(x-x_{j})+(\bm{W}_{y})_{j,k}(y-y_{k}),\quad(x,y)\in C_{j,k},

(2.10)

which results in

		$\displaystyle\bm{W}^{\rm E}_{j,k}=\bm{W}_{j,k}+(\bm{W}_{x})_{j,k}\frac{\Delta x}{2},\quad\bm{W}^{\rm W}_{j,k}=\bm{W}_{j,k}-(\bm{W}_{x})_{j,k}\frac{\Delta x}{2},$		(2.11)
		$\displaystyle\bm{W}^{\rm N}_{j,k}=\bm{W}_{j,k}+(\bm{W}_{y})_{j,k}\frac{\Delta y}{2},\quad\bm{W}^{\rm S}_{j,k}=\bm{W}_{j,k}-(\bm{W}_{y})_{j,k}\frac{\Delta y}{2}.$		(2.11)

In order for (2.10) to be second-order accurate, the slopes $(\bm{W}_{x})_{j,k}$ and $(\bm{W}_{y})_{j,k}$ have to be at least first-order approximations of $\bm{W}_{x}(x_{j},y_{k})$ and $\bm{W}_{y}(x_{j},y_{k})$ , respectively. A non-oscillatory nature of the piecewise linear reconstruction (2.10) is typically ensured with the help of a nonlinear limiter. To all of the components of $\bm{W}$ , we compute the slopes (except for $((b_{1})_{x})_{j,k}$ and $((b_{2})_{y})_{j,k}$ ) using the generalized minmod limiter (see, e.g., [41, 46, 55]):

		$\displaystyle(W^{(i)}_{x})_{j,k}={\rm minmod}\left(\theta\,\frac{W^{(i)}_{j,k}-W^{(i)}_{j-1,k}}{\Delta x},\,\frac{W^{(i)}_{j+1,k}-W^{(i)}_{j-1,k}}{2\Delta x},\,\theta\,\frac{W^{(i)}_{j+1,k}-W^{(i)}_{j,k}}{\Delta x}\right),\quad i\neq 5,$		(2.12)
		$\displaystyle(W^{(i)}_{y})_{j,k}={\rm minmod}\left(\theta\,\frac{W^{(i)}_{j,k}-W^{(i)}_{j,k-1}}{\Delta y},\,\frac{W^{(i)}_{j,k+1}-W^{(i)}_{j,k-1}}{2\Delta y},\,\theta\,\frac{W^{(i)}_{j,k+1}-W^{(i)}_{j,k}}{\Delta y}\right),\quad i\neq 6,$		(2.12)

where the minmod function is defined by

{\rm minmod}(a_{1},a_{2},\dots)=\begin{cases}{\rm sgn}(a_{k})\min\{{\left|a_{1}\right|},{\left|a_{2}\right|},\dots\}&\mbox{if}\leavevmode\nobreak\ {\rm sgn}(a_{k})=1\leavevmode\nobreak\ \mbox{or}\leavevmode\nobreak\ {\rm sgn}(a_{k})=-1\leavevmode\nobreak\ \forall k,\\ 0&\mbox{otherwise}.\end{cases}

(2.13)

The slopes $((b_{1})_{x})_{j,k}$ and $((b_{2})_{y})_{j,k}$ , however, should not be computed using the generalized minmod limiter or any other conventional limiter as our goal is to enforce local divergence-free condition (2.3), which at the discrete level reads as $((b_{1})_{x})_{j,k}+((b_{2})_{y})_{j,k}\equiv 0$ for all $j,k$ . This goal can be achieved if we set

((b_{1})_{x})_{j,k}=\,\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$A$\kern 0.0pt}}}_{j,k}\quad\mbox{and}\quad((b_{1})_{y})_{j,k}=\,\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$B$\kern 0.0pt}}}_{j,k},

(2.14)

since

\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$A$\kern 0.0pt}}}_{j,k}+\,\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$B$\kern 0.0pt}}}_{j,k}=0

(2.15)

is true for all $j,k$ provided (2.15) is satisfied at time $t=0$ ; see Theorem 2.2 in §2.2.4.

While the use of (2.14) guarantees the local discrete divergence-free condition, the resulting reconstruction of $b_{1}$ and $b_{2}$ may be oscillatory in the $x$ - and $y$ -directions, respectively. As we have observed in several numerical experiments, this often leads to an oscillatory numerical solution. We, therefore, adjust the slopes in (2.14) by scaling them as follows.

We begin by introducing the auxiliary slopes $(\widehat{(b_{1})}_{x})_{j,k}$ and $(\widehat{(b_{2})}_{y})_{j,k}$ , which are computed using the aforementioned generalized minmod reconstruction. The reconstructions of $b_{1}$ and $b_{2}$ can then be made both non-oscillatory and locally divergence-free by replacing (2.14) with

((b_{1})_{x})_{j,k}=\sigma_{j,k}\,\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$A$\kern 0.0pt}}}_{j,k},\quad((b_{1})_{y})_{j,k}=\sigma_{j,k}\,\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$B$\kern 0.0pt}}}_{j,k},

(2.16)

where

\sigma_{j,k}=\min\Big{\{}1,\sigma^{x}_{j,k},\sigma^{y}_{j,k}\Big{\}}

(2.17)

and the scaling factors $\sigma^{x}_{j,k}$ and $\sigma^{y}_{j,k}$ are computed by

\sigma^{x}_{j,k}:=\left\{\begin{aligned} &\min\bigg{\{}1,\frac{(\widehat{(b_{1})}_{x})_{j,k}}{\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$A$\kern 0.0pt}}}_{j,k}}\bigg{\}}&&\mbox{if}\leavevmode\nobreak\ (\widehat{(b_{1})}_{x})_{j,k}\,\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$A$\kern 0.0pt}}}_{j,k}>0,\\ &0&&\mbox{otherwise},\end{aligned}\right.

(2.18)

and

\sigma^{y}_{j,k}:=\left\{\begin{aligned} &\min\bigg{\{}1,\frac{(\widehat{(b_{2})}_{y})_{j,k}}{\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$B$\kern 0.0pt}}}_{j,k}}\bigg{\}}&&\mbox{if}\leavevmode\nobreak\ (\widehat{(b_{2})}_{y})_{j,k}\,\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$B$\kern 0.0pt}}}_{j,k}>0,\\ &0&&\mbox{otherwise}.\end{aligned}\right.

(2.19)

Finally, equipped with (2.11), we use the EOS (2.2) to compute the cell interface point values ${\cal E}^{\rm E,\rm W,\rm N,\rm S}_{j,k}$ as follows:

{\cal E}_{j,k}^{\ell}=\frac{p_{j,k}^{\ell}}{\gamma-1}+\frac{1}{2}\rho^{\ell}_{j,k}\left[(u^{\ell}_{j,k})^{2}+(v^{\ell}_{j,k})^{2}+(w^{\ell}_{j,k})^{2}\right]+\frac{1}{2}\left[{\left((b_{1})^{\ell}_{j,k}\right)}^{2}+{\left((b_{2})^{\ell}_{j,k}\right)}^{2}+{\left((b_{3})^{\ell}_{j,k}\right)}^{2}\right],

(2.20)

where $\ell\in\{{\rm E},{\rm W},{\rm N},{\rm S}\}$ .

Remark 2.1.

We note that we have reconstructed the primitive variables $u$ , $v$ , $w$ , and $p$ rather than the conservative variables $\rho u$ , $\rho v$ , $\rho w$ , and ${\cal E}$ since our numerical experiments clearly indicate that the resulting scheme, which is based on the reconstruction of the primitive variables, is less oscillatory and produces no negative pressure values.

It is important to point out that the ninth and tenth components of the fluxes $\bm{F}(\bm{U}^{\rm E(W)}_{j,k})$ and $\bm{G}(\bm{U}^{\rm N(S)}_{j,k})$ depend not only on the corresponding point values of $u$ , $v$ , $A$ , $B$ , $b_{1}$ , and $b_{2}$ , but also on the point values of the derivatives $(u_{y})^{\rm E(W)}_{j,k}$ and $(v_{x})^{\rm N(S)}_{j,k}$ . We compute these values using the first-order approximation, namely, we set

(u_{y})^{\rm E}_{j,k}=(u_{y})_{j,k},\quad(u_{y})^{\rm W}_{j,k}=(u_{y})_{j,k},\quad(v_{x})^{\rm N}_{j,k}=(v_{x})_{j,k},\quad(v_{x})^{\rm S}_{j,k}=(v_{x})_{j,k},

where the slopes $(u_{y})_{j,k}$ and $(v_{x})_{j,k}$ are computed by (2.12). Notice that even though this will result in the first-order approximation of the auxiliary variables $A$ and $B$ , the other components of $\bm{U}$ will be still computed with the second order and thus the second-order accuracy of the resulting scheme will not be affected.

2.2.2 One-Sided Speeds of Propagation

Equipped with the reconstructed point values (2.11) and (2.20), we now proceed with the computation of the one-sided local speeds of propagation $s^{\pm}_{{j+\frac{1}{2}},k}$ and $s^{\pm}_{j,{k+\frac{1}{2}}}$ seen in (2.8) and (2.9). We stress that when the PCCU schemes are applied to general nonconservative systems of type (2.6), the $x$ - and $y$ -directional speeds would typically be estimated using the largest and smallest eigenvalues of the matrices $\frac{\partial\bm{F}}{\partial\bm{U}}(\bm{U})-Q(\bm{U})$ and $\frac{\partial\bm{G}}{\partial\bm{U}}(\bm{U})-R(\bm{U})$ , respectively. However, it is known (see, e.g., [20]) that in the context of the ideal MHD system (2.6)–(LABEL:2.10), the estimates, which are solely based on the eigenvalues mentioned above may be inaccurate and using them may lead to severe instabilities.

We, therefore, follow [63], where the propagation speeds were slightly overestimated to ensure the positivity of both the computed density and pressure, and estimate the right- and left-sided local speeds in the $x$ -direction by

		$\displaystyle s_{{j+\frac{1}{2}},k}^{+}=\max\Big{\{}\max\big{\{}u_{j,k}^{\rm E},u_{{j+\frac{1}{2}},k}^{\rm Roe}\big{\}}+c_{j,k}^{\rm E}+\beta^{x}_{{j+\frac{1}{2}},k},\,\max\big{\{}u_{j+1,k}^{\rm W},u_{{j+\frac{1}{2}},k}^{\rm Roe}\big{\}}+c_{j+1,k}^{\rm W}+\beta^{x}_{{j+\frac{1}{2}},k},\,0\Big{\}},$
		$\displaystyle s_{{j+\frac{1}{2}},k}^{-}=\min\Big{\{}\min\big{\{}u_{j,k}^{\rm E},u_{{j+\frac{1}{2}},k}^{\rm Roe}\big{\}}-c_{j,k}^{\rm E}-\beta^{x}_{{j+\frac{1}{2}},k},\,\min\big{\{}u_{j+1,k}^{\rm W},u_{{j+\frac{1}{2}},k}^{\rm Roe}\big{\}}-c_{j+1,k}^{\rm W}-\beta^{x}_{{j+\frac{1}{2}},k},\,0\Big{\}},$

where

u_{{j+\frac{1}{2}},k}^{\rm Roe}:=\frac{u_{j,k}^{\rm E}\sqrt{\rho_{j,k}^{\rm E}}+u_{j+1,k}^{\rm W}\sqrt{\rho_{j+1,k}^{\rm W}}}{\sqrt{\rho_{j,k}^{\rm E}}+\sqrt{\rho_{j+1,k}^{\rm W}}},\quad\beta^{x}_{{j+\frac{1}{2}},k}:=\frac{{\left|\bm{b}_{j,k}^{\rm E}-\bm{b}_{j+1,k}^{\rm W}\right|}}{\sqrt{\rho_{j,k}^{\rm E}}+\sqrt{\rho_{j+1,k}^{\rm W}}},

and $c^{\rm E(W)}_{j,k}$ are the fast magneto-acoustic wave speeds computed using

\left(c_{j,k}^{\rm E(W)}\right)^{2}=\frac{1}{2\rho_{j,k}^{\rm E(W)}}{\left[\gamma p_{j,k}^{\rm E(W)}+{\left|\bm{b}_{j,k}^{\rm E(W)}\right|}^{2}+\sqrt{{\left(\gamma p_{j,k}^{\rm E(W)}+{\left|\bm{b}_{j,k}^{\rm E(W)}\right|}^{2}\right)}^{2}-4\gamma p_{j,k}^{\rm E(W)}{\left({\left(b_{1}\right)}_{j,k}^{\rm E(W)}\right)}^{2}}\leavevmode\nobreak\ \right]}.

Similarly, we estimate the corresponding $y$ -directional speeds by

		$\displaystyle s_{j,{k+\frac{1}{2}}}^{+}=\max\Big{\{}\max\big{\{}v_{j,k}^{\rm N},v_{j,{k+\frac{1}{2}}}^{\rm Roe}\big{\}}+c_{j,k}^{\rm N}+\beta^{y}_{j,{k+\frac{1}{2}}},\,\max\big{\{}v_{j,k+1}^{\rm S},v_{j,{k+\frac{1}{2}}}^{\rm Roe}\big{\}}+c_{j,k+1}^{\rm S}+\beta^{y}_{j,{k+\frac{1}{2}}},\,0\Big{\}},$
		$\displaystyle s_{j,{k+\frac{1}{2}}}^{-}=\min\Big{\{}\min\big{\{}v_{j,k}^{\rm N},v_{j,{k+\frac{1}{2}}}^{\rm Roe}\big{\}}-c_{j,k}^{\rm N}-\beta^{y}_{j,{k+\frac{1}{2}}},\,\min\big{\{}v_{j,k+1}^{\rm S},v_{j,{k+\frac{1}{2}}}^{\rm Roe}\big{\}}-c_{j,k+1}^{\rm S}-\beta^{y}_{j,{k+\frac{1}{2}}},\,0\Big{\}},$

where

v_{j,{k+\frac{1}{2}}}^{\rm Roe}:=\frac{v_{j,k}^{\rm N}\sqrt{\rho_{j,k}^{\rm N}}+v_{j,k+1}^{\rm S}\sqrt{\rho_{j,k+1}^{\rm S}}}{\sqrt{\rho_{j,k}^{\rm N}}+\sqrt{\rho_{j,k+1}^{\rm S}}},\quad\beta^{y}_{j,{k+\frac{1}{2}}}:=\frac{{\left|\bm{b}_{j,k}^{\rm N}-\bm{b}_{j,k+1}^{\rm S}\right|}}{\sqrt{\rho_{j,k}^{\rm N}}+\sqrt{\rho_{j,k+1}^{\rm S}}},

\left(c_{j,k}^{\rm N(S)}\right)^{2}=\frac{1}{2\rho_{j,k}^{\rm N(S)}}{\left[\gamma p_{j,k}^{\rm N(S)}+{\left|\bm{b}_{j,k}^{\rm N(S)}\right|}^{2}+\sqrt{{\left(\gamma p_{j,k}^{\rm N(S)}+{\left|\bm{b}_{j,k}^{\rm N(S)}\right|}^{2}\right)}^{2}-4\gamma p_{j,k}^{\rm N(S)}{\left({\left(b_{1}\right)}_{j,k}^{\rm N(S)}\right)}^{2}}\leavevmode\nobreak\ \right]}.

2.2.3 Discretization of the Nonconservative Products

In this section, we provide the computation of the nonconservative product terms in (2.8).

Following [8] (see also [11]), we obtain nonconservative terms in the $x$ -direction, $\bm{Q}_{j,k}$ and $\bm{Q}_{\bm{\Psi},{j+\frac{1}{2}},k}$ , as follows. First, in order to compute the term $\bm{Q}_{j,k}$ , we take a global (in space) interpolant $\bm{U}(\widetilde{\bm{W}}(x,y))$ , where $\widetilde{\bm{W}}$ is given by (2.10), and evaluate the integral in

\bm{Q}_{j,k}=\int\limits_{x_{j-\frac{1}{2}}}^{x_{j+\frac{1}{2}}}Q\big{(}\bm{U}\big{(}\widetilde{\bm{W}}(x,y_{k})\big{)}\big{)}\bm{U}(\widetilde{\bm{W}}(x,y_{k}))_{x}\,{\rm d}x,

where $Q(\bm{U})$ is defined in (LABEL:2.10), exactly. This results in the following expressions for the ten components of the vector $\bm{Q}_{j,k}$ :

	$\displaystyle Q^{(1)}_{j,k}=Q^{(9)}_{j,k}=Q^{(10)}_{j,k}=0,$
	$\displaystyle Q^{(i)}_{j,k}=-\int\limits_{x_{j-\frac{1}{2}}}^{x_{j+\frac{1}{2}}}\widetilde{b_{i-1}}(x,y_{k})((b_{1})_{x})_{j,k}\,{\rm d}x=-(\,\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$b_{i-1}$\kern 0.0pt}}})_{j,k}\sigma_{j,k}\,\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$A$\kern 0.0pt}}}_{j,k}\Delta x,\quad i=2,3,4,$
	$\displaystyle\big{(}Q^{(5)}_{j,k},Q^{(6)}_{j,k},Q^{(7)}_{j,k}\big{)}^{\top}=-\int\limits_{x_{j-\frac{1}{2}}}^{x_{j+\frac{1}{2}}}\widetilde{\bm{u}}(x,y_{k})((b_{1})_{x})_{j,k}\,{\rm d}x=-\bm{u}_{j,k}\sigma_{j,k}\,\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$A$\kern 0.0pt}}}_{j,k}\Delta x,$
	$\displaystyle Q^{(8)}_{j,k}=-\int\limits_{x_{j-\frac{1}{2}}}^{x_{j+\frac{1}{2}}}\widetilde{\bm{u}}(x,y_{k})\cdot\widetilde{\bm{b}}(x,y_{k})((b_{1})_{x})_{j,k}\,{\rm d}x=-\Big{[}\Big{(}\bm{u}_{j,k}\cdot\bm{b}_{j,k}$
	$\displaystyle\hskip 51.21504pt+\frac{(\Delta x)^{2}}{12}\left\{(u_{x})_{j,k}\sigma_{j,k}\,\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$A$\kern 0.0pt}}}_{j,k}+(v_{x})_{j,k}((b_{2})_{x})_{j,k}+(w_{x})_{j,k}((b_{3})_{x})_{j,k}\right\}\Big{)}\Big{]}\sigma_{j,k}\,\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$A$\kern 0.0pt}}}_{j,k}\Delta x,$

where we have used the slopes $((b_{1})_{x})_{j,k}$ given by (2.16)–(2.19), while the other slopes are computed in (2.12)–(2.13).

Next, the terms $\bm{Q}_{\bm{\Psi},{j+\frac{1}{2}},k}$ are computed by the exact integration of

\bm{Q}_{\bm{\Psi},{j+\frac{1}{2}},k}=\int\limits_{0}^{1}Q\big{(}\bm{U}\big{(}\bm{\Psi}_{{j+\frac{1}{2}},k}(s)\big{)}\big{)}\bm{\Psi}_{{j+\frac{1}{2}},k}^{\prime}(s)\,{\rm d}s,

where $\bm{\Psi}_{{j+\frac{1}{2}},k}(s)$ is a linear path connecting the states $\bm{W}^{\rm E}_{j,k}$ and $\bm{W}^{\rm W}_{j+1,k}$ :

\bm{\Psi}_{{j+\frac{1}{2}},k}(s)=\bm{W}^{\rm E}_{j,k}+s{\left(\bm{W}^{\rm W}_{j+1,k}-\bm{W}^{\rm E}_{j,k}\right)}.

This results in

	$\displaystyle Q^{(1)}_{\bm{\Psi},{j+\frac{1}{2}},k}=Q^{(9)}_{\bm{\Psi},{j+\frac{1}{2}},k}=Q^{(10)}_{\bm{\Psi},{j+\frac{1}{2}},k}=0,$
	$\displaystyle Q^{(i)}_{\bm{\Psi},{j+\frac{1}{2}},k}=-\int\limits_{0}^{1}\left\{(b_{i-1})^{\rm E}_{j,k}+s{\left((b_{i-1})^{\rm W}_{j+1,k}-(b_{i-1})^{\rm E}_{j,k}\right)}\right\}[b_{1}]_{{j+\frac{1}{2}},k}\,{\rm d}s$
	$\displaystyle\hskip 44.10185pt=-\frac{1}{2}\left((b_{i-1})^{\rm E}_{j,k}+(b_{i-1})^{\rm W}_{j+1,k}\right)[b_{1}]_{{j+\frac{1}{2}},k},\quad i=2,3,4,$
	$\displaystyle\big{(}Q^{(5)}_{\bm{\Psi},{j+\frac{1}{2}},k},Q^{(6)}_{\bm{\Psi},{j+\frac{1}{2}},k},Q^{(7)}_{\bm{\Psi},{j+\frac{1}{2}},k}\big{)}^{\top}=-\int\limits_{0}^{1}\left\{\bm{u}^{\rm E}_{j,k}+s{\left(\bm{u}^{\rm W}_{j+1,k}-\bm{u}^{\rm E}_{j,k}\right)}\right\}[b_{1}]_{{j+\frac{1}{2}},k}\,{\rm d}s$
	$\displaystyle\hskip 159.3356pt=-\frac{1}{2}\left(\bm{u}^{\rm E}_{j,k}+\bm{u}^{\rm W}_{j+1,k}\right)[b_{1}]_{{j+\frac{1}{2}},k},$
	$\displaystyle Q^{(8)}_{\bm{\Psi},{j+\frac{1}{2}},k}=-\int\limits_{0}^{1}\left\{\bm{u}^{\rm E}_{j,k}+s{\left(\bm{u}^{\rm W}_{j+1,k}-\bm{u}^{\rm E}_{j,k}\right)}\right\}\cdot\left\{\bm{b}^{\rm E}_{j,k}+s{\left(\bm{b}^{\rm W}_{j+1,k}-\bm{b}^{\rm E}_{j,k}\right)}\right\}[b_{1}]_{{j+\frac{1}{2}},k}\,{\rm d}s$
	$\displaystyle\hskip 44.10185pt=-\frac{1}{6}\Big{(}2\bm{u}_{j,k}^{\rm E}\cdot\bm{b}_{j,k}^{\rm E}+\bm{u}_{j,k}^{\rm E}\cdot\bm{b}_{j+1,k}^{\rm W}+\bm{u}_{j+1,k}^{\rm W}\cdot\bm{b}_{j,k}^{\rm E}+2\bm{u}_{j+1,k}^{\rm W}\cdot\bm{b}_{j+1,k}^{\rm W}\Big{)}[b_{1}]_{{j+\frac{1}{2}},k},$

where $[b_{1}]_{{j+\frac{1}{2}},k}:=(b_{1})^{\rm W}_{j+1,k}-(b_{1})^{\rm E}_{j,k}$ .

Similarly, we obtain the following formulae for the nonconservative terms in the $y$ -direction, $\bm{R}_{j,k}$ and $\bm{R}_{\bm{\Psi},j,{k+\frac{1}{2}}}$ :

	$\displaystyle R^{(1)}_{j,k}=R^{(9)}_{j,k}=R^{(10)}_{j,k}=0,\quad R^{(i)}_{j,k}=-(\,\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$b_{i-1}$\kern 0.0pt}}})_{j,k}\sigma_{j,k}\,\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$B$\kern 0.0pt}}}_{j,k}\Delta y,\quad i=2,3,4,$
	$\displaystyle\big{(}R^{(5)}_{j,k},R^{(6)}_{j,k},R^{(7)}_{j,k}\big{)}^{\top}=-\bm{u}_{j,k}\sigma_{j,k}\,\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$B$\kern 0.0pt}}}_{j,k}\Delta y,$
	$\displaystyle R^{(8)}_{j,k}=-\Big{[}\Big{(}\bm{u}_{j,k}\cdot\bm{b}_{j,k}+\frac{(\Delta y)^{2}}{12}\left\{(u_{y})_{j,k}((b_{1})_{y})_{j,k}+(v_{y})_{j,k}\sigma_{j,k}\,\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$B$\kern 0.0pt}}}_{j,k}+(w_{y})_{j,k}((b_{3})_{y})_{j,k}\right\}\Big{)}\Big{]}\sigma_{j,k}\,\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$B$\kern 0.0pt}}}_{j,k}\Delta y,$

where we have used the slopes $((b_{2})_{y})_{j,k}$ given by (2.16)–(2.19), while the other slopes are computed in (2.12)–(2.13), and

	$\displaystyle R^{(1)}_{\bm{\Psi},j,{k+\frac{1}{2}}}=R^{(9)}_{\bm{\Psi},j,{k+\frac{1}{2}}}=R^{(10)}_{\bm{\Psi},j,{k+\frac{1}{2}}}=0,\quad R^{(i)}_{\bm{\Psi},j,{k+\frac{1}{2}}}=-\frac{1}{2}\left((b_{i-1})^{\rm N}_{j,k}+(b_{i-1})^{\rm S}_{j,k+1}\right)[b_{2}]_{j,{k+\frac{1}{2}}},\quad i=2,3,4,$
	$\displaystyle\big{(}R^{(5)}_{\bm{\Psi},j,{k+\frac{1}{2}}},R^{(6)}_{\bm{\Psi},j,{k+\frac{1}{2}}},R^{(7)}_{\bm{\Psi},j,{k+\frac{1}{2}}}\big{)}^{\top}=-\frac{1}{2}\left(\bm{u}^{\rm N}_{j,k}+\bm{u}^{\rm S}_{j,k+1}\right)[b_{2}]_{j,{k+\frac{1}{2}}},$
	$\displaystyle R^{(8)}_{\bm{\Psi},j,{k+\frac{1}{2}}}=-\frac{1}{6}\Big{(}2\bm{u}_{j,k}^{\rm N}\cdot\bm{b}_{j,k}^{\rm N}+\bm{u}_{j,k}^{\rm N}\cdot\bm{b}_{j,k+1}^{\rm S}+\bm{u}_{j,k+1}^{\rm S}\cdot\bm{b}_{j,k}^{\rm N}+2\bm{u}_{j,k+1}^{\rm S}\cdot\bm{b}_{j,k+1}^{\rm S}\Big{)}[b_{2}]_{j,{k+\frac{1}{2}}},$

where $[b_{2}]_{j,{k+\frac{1}{2}}}:={\left((b_{2})^{\rm S}_{j,k+1}-(b_{2})^{\rm N}_{j,k}\right)}$ .

2.2.4 Local Divergence-Free Property

We now prove the local divergence-free property of the proposed PCCU scheme.

Theorem 2.2.

For the PCCU scheme (2.8)–(2.9) with the reconstruction described in §2.2.1, the local divergence-free condition

\big{(}(b_{1})_{x}\big{)}_{j,k}+\big{(}(b_{2})_{y}\big{)}_{j,k}=0,

(2.21)

holds for all $j,k$ and at all times, provided it is satisfied initially.

Proof.

First, we note that according to (2.16),

\big{(}(b_{1})_{x}\big{)}_{j,k}+\big{(}(b_{2})_{y}\big{)}_{j,k}=\sigma_{j,k}\left(\,\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$A$\kern 0.0pt}}}_{j,k}+\,\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$B$\kern 0.0pt}}}_{j,k}\right).

Therefore, in order to prove (2.21), it is sufficient to show that $\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$A$\kern 0.0pt}}}_{j,k}+\,\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$B$\kern 0.0pt}}}_{j,k}=0$ for all $j,k$ and for all times assuming that it is satisfied at the initial time $t=0$ .

We then observe that the quantities $\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$A$\kern 0.0pt}}}_{j,k}$ and $\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$B$\kern 0.0pt}}}_{j,k}$ are the ninth and tenth components of $\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$\bm{U}$\kern 0.0pt}}}_{j,k}$ and thus they are evolved in time by numerically integrating the ninth and tenth components of (2.8)–(2.9). Adding these components in (2.8) results in

	$\displaystyle\frac{{\rm d}}{{\rm d}t}\left(\,\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$A$\kern 0.0pt}}}_{j,k}+\,\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$B$\kern 0.0pt}}}_{j,k}\right)=$	$\displaystyle-\frac{1}{\Delta x}\bigg{[}{{\cal F}}^{(9)}_{{j+\frac{1}{2}},k}-{{\cal F}}^{(9)}_{{j-\frac{1}{2}},k}+{{\cal F}}^{(10)}_{{j+\frac{1}{2}},k}-{{\cal F}}^{(10)}_{{j-\frac{1}{2}},k}\bigg{]}$		(2.22)
		$\displaystyle-\frac{1}{\Delta y}\bigg{[}{{\cal G}}^{(9)}_{j,{k+\frac{1}{2}}}-{{\cal G}}^{(9)}_{j,{k-\frac{1}{2}}}+{{\cal G}}^{(10)}_{j,{k+\frac{1}{2}}}-{{\cal G}}^{(10)}_{j,{k-\frac{1}{2}}}\bigg{]}.$		(2.22)

In order to complete the proof, it is sufficient to show that the RHS of (2.22) vanishes as long as $\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$A$\kern 0.0pt}}}_{j,k}+\,\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$B$\kern 0.0pt}}}_{j,k}=0$ for all $j,k$ . To this end, we use (2.9) to evaluate

	$\displaystyle{{\cal F}}^{(9)}_{{j+\frac{1}{2}},k}-{{\cal F}}^{(9)}_{{j-\frac{1}{2}},k}+{{\cal F}}^{(10)}_{{j+\frac{1}{2}},k}-{{\cal F}}^{(10)}_{{j-\frac{1}{2}},k}$	$\displaystyle=\frac{s_{{j+\frac{1}{2}},k}^{+}\left[u^{\rm E}_{j,k}\left(A_{j,k}^{\rm E}+B_{j,k}^{\rm E}\right)\right]-s_{{j+\frac{1}{2}},k}^{-}\left[u^{\rm W}_{j+1,k}\left(A_{j+1,k}^{\rm W}+B_{j+1,k}^{\rm W}\right)\right]}{s_{{j+\frac{1}{2}},k}^{+}-s_{{j+\frac{1}{2}},k}^{-}}$
		$\displaystyle+\frac{s_{{j+\frac{1}{2}},k}^{+}s_{{j+\frac{1}{2}},k}^{-}}{s_{{j+\frac{1}{2}},k}^{+}-s_{{j+\frac{1}{2}},k}^{-}}\Big{[}\left(A^{\rm W}_{j+1,k}+B^{\rm W}_{j+1,k}\right)-\left(A^{\rm E}_{j,k}+B^{\rm E}_{j,k}\right)\Big{]}$
	$\displaystyle\stackrel{{\scriptstyle(\ref{2.15})}}{{=}}\frac{s_{{j+\frac{1}{2}},k}^{+}}{s_{{j+\frac{1}{2}},k}^{+}-s_{{j+\frac{1}{2}},k}^{-}}$	$\displaystyle\left[u^{\rm E}_{j,k}\Big{(}\,\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$A$\kern 0.0pt}}}_{j,k}+\frac{\Delta x}{2}(A_{x})_{j,k}+\,\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$B$\kern 0.0pt}}}_{j,k}+\frac{\Delta x}{2}(B_{x})_{j,k}\Big{)}\right]$
	$\displaystyle-\frac{s_{{j+\frac{1}{2}},k}^{-}}{s_{{j+\frac{1}{2}},k}^{+}-s_{{j+\frac{1}{2}},k}^{-}}$	$\displaystyle\left[u^{\rm W}_{j+1,k}\Big{(}\,\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$A$\kern 0.0pt}}}_{j+1,k}-\frac{\Delta x}{2}(A_{x})_{j+1,k}+\,\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$B$\kern 0.0pt}}}_{j+1,k}-\frac{\Delta x}{2}(B_{x})_{j+1,k}\Big{)}\right]$
	$\displaystyle+\frac{s_{{j+\frac{1}{2}},k}^{+}s_{{j+\frac{1}{2}},k}^{-}}{s_{{j+\frac{1}{2}},k}^{+}-s_{{j+\frac{1}{2}},k}^{-}}$	$\displaystyle\left[\,\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$A$\kern 0.0pt}}}_{j+1,k}-\frac{\Delta x}{2}(A_{x})_{j+1,k}+\,\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$B$\kern 0.0pt}}}_{j+1,k}-\frac{\Delta x}{2}(B_{x})_{j+1,k}\right]$
	$\displaystyle-\frac{s_{{j+\frac{1}{2}},k}^{+}s_{{j+\frac{1}{2}},k}^{-}}{s_{{j+\frac{1}{2}},k}^{+}-s_{{j+\frac{1}{2}},k}^{-}}$	$\displaystyle\left[\,\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$A$\kern 0.0pt}}}_{j,k}+\frac{\Delta x}{2}(A_{x})_{j,k}+\,\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$B$\kern 0.0pt}}}_{j,k}+\frac{\Delta x}{2}(B_{x})_{j,k}\right]$
		$\displaystyle=\frac{\Delta x}{2}\Bigg{\{}\frac{s_{{j+\frac{1}{2}},k}^{+}\left[u^{\rm E}_{j,k}\left((A_{x})_{j,k}+(B_{x})_{j,k}\right)\right]+s_{{j+\frac{1}{2}},k}^{-}\left[u^{\rm W}_{j+1,k}\left((A_{x})_{j+1,k}+(B_{x})_{j+1,k}\right)\right]}{s_{{j+\frac{1}{2}},k}^{+}-s_{{j+\frac{1}{2}},k}^{-}}$
		$\displaystyle-\frac{s_{{j+\frac{1}{2}},k}^{+}s_{{j+\frac{1}{2}},k}^{-}}{s_{{j+\frac{1}{2}},k}^{+}-s_{{j+\frac{1}{2}},k}^{-}}\Big{[}(A_{x})_{j+1,k}+(B_{x})_{j+1,k}+(A_{x})_{j,k}+(B_{x})_{j,k}\Big{]}\Bigg{\}},$

where the last equality is obtained using $\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$A$\kern 0.0pt}}}_{j,k}+\,\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$B$\kern 0.0pt}}}_{j,k}=\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$A$\kern 0.0pt}}}_{j+1,k}+\,\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$B$\kern 0.0pt}}}_{j+1,k}=0$ .

It is now clear that ${{\cal F}}^{(9)}_{{j+\frac{1}{2}},k}-{{\cal F}}^{(9)}_{{j-\frac{1}{2}},k}+{{\cal F}}^{(10)}_{{j+\frac{1}{2}},k}-{{\cal F}}^{(10)}_{{j-\frac{1}{2}},k}$ will be identically zero as long as

(A_{x})_{j,k}+(B_{x})_{j,k}=0

(2.23)

for all $j,k$ . Indeed, (2.23) is true since $\,\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$A$\kern 0.0pt}}}_{j,k}=-\,\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$B$\kern 0.0pt}}}_{j,k}$ for all $j,k$ and both the slopes $(A_{x})_{j,k}$ and $(B_{x})_{j,k}$ are computed using the same limiter (2.12).

Similarly, one can show that ${{\cal G}}^{(9)}_{j,{k+\frac{1}{2}}}-{{\cal G}}^{(9)}_{j,{k-\frac{1}{2}}}+{{\cal G}}^{(10)}_{j,{k+\frac{1}{2}}}-{{\cal G}}^{(10)}_{j,{k-\frac{1}{2}}}\equiv 0$ , so that the RHS of (2.22) vanishes and thus the proof of the theorem is complete. ∎

2.3 Numerical Examples

In this section, we demonstrate the performance of the proposed PCCU scheme in several numerical experiments conducted for the augmented 2-D ideal MHD system (2.6)-(LABEL:2.10), (2.2). In all of the examples in this section, we take the CFL number 0.25 and the minmod parameter $\theta=1.3$ .

Example 1—Brio-Wu Shock-Tube Problem.

In the first example, we consider the one-dimensional (1-D) Riemann problem known as the Brio-Wu shock tube problem, originally presented in [7]. This problem is the standard test for capturing compound waves that emerge as solutions of the ideal MHD system. We take the following initial data, which depend on $x$ only:

\left(\rho,u,v,w,b_{1},b_{2},b_{3},p\right)(x,y,0)=\left\{\begin{aligned} &(1,0,0,0,0.75,1,0,1),&&x<0,\\ &(0.125,0,0,0,0.75,-1,0,0.1),&&\mbox{otherwise},\end{aligned}\right.

and set the free boundary conditions on all sides of the computational domain $[-1,1]\times[-0.01,0.01]$ . The specific heat ratio $\gamma=2$ in this example.

In Figure 2.1, we plot the $y=0$ cross-section of the density $\rho$ , $x$ -magnetic field $b_{1}$ , and $y$ -magnetic field $b_{2}$ computed on $800\times 8$ and $1600\times 16$ uniform meshes at time $t=0.2$ . As one can see, the solution to this 1-D Riemann problem consists of several nonsmooth structures, such as rarefaction waves, shocks traveling at various speeds, a contact discontinuity, and a compound shock wave. The proposed PCCU scheme captures all of these complicated structures well, and the obtained results strongly agree with those reported in [19, 23, 38, 40, 43].

Refer to caption — Figure 2.1: Example 1: $\rho$ , $b_{1}$ , and $b_{2}$ computed by the PCCU scheme on $800\times 8$ (red circles) and $1600\times 16$ (solid black line) uniform meshes.

Example 2—Orszag-Tang Vortex Problem.

In the second example, we consider the Orszag-Tang vortex problem, which was introduced in [47] and has been widely used as a benchmark due to the formation and interaction of multiple shocks as the system evolves in time and the presence of many important features of MHD turbulence; see, e.g., [19, 40, 39, 42, 65]. The initial conditions for this problem read as

		$\displaystyle\rho(x,y,0)\equiv\gamma^{2},\quad u(x,y,0)=-\sin y,\quad v(x,y,0)=\sin x,\quad w(x,y,0)\equiv 0,$
		$\displaystyle b_{1}(x,y,0)=-\sin y,\quad b_{2}(x,y,0)=\sin(2x),\quad b_{3}(x,y,0)\equiv 0,\quad p(x,y,0)\equiv\gamma,$

where $\gamma=5/3$ is the specific heat ratio. We set the periodic boundary conditions on all sides of the computational domain $[0,2\pi]\times[0,2\pi]$ .

The time evolution of the fluid density $\rho$ computed on a uniform $200\times 200$ mesh is shown at times $t=0.5$ , 2, 3, and 4 in Figure 2.2. We observe that the numerical solution computed by the proposed PCCU scheme remains stable and is consistent with previous results presented in [40, 39, 42, 65], demonstrating that the ability of our scheme to capture both smooth flows and shocks.

Example 3—Rotor Problem.

Next, we consider the “second rotor problem” from [5, 56], referred to as the rotor problem in this paper. This commonly used benchmark problem describes a rapidly-rotating disk of dense fluid centered in a background of static fluid. Over time, the disk expands and rotates. The initial conditions are given by

		$\displaystyle(\rho,u,v)(x,y,0)=\left\{\begin{aligned} &\left(10,\frac{0.5-y}{r_{0}},\frac{x-0.5}{r_{0}}\right),&&r<0.1,\\ &\left(1+9\lambda,\frac{\lambda(0.5-y)}{r},\frac{\lambda(x-0.5)}{r}\right),&&0.1\leq r\leq 0.115,\\ &(1,0,0),&&r>0.115,\end{aligned}\right.$
		$\displaystyle w(x,y,0)=b_{2}(x,y,0)=b_{3}(x,y,0)\equiv 0,\quad b_{1}(x,y,0)\equiv\frac{2.5}{\sqrt{4\pi}},\quad p(x,y,0)\equiv 0.5,$

where $r=\sqrt{(x-0.5)^{2}+(y-0.5)^{2}}$ and $\lambda=(0.115-r)/0.015$ . We take the specific heat ratio $\gamma=5/3$ and use the periodic boundary conditions on all sides of the computational domain $[0,1]\times[0,1]$ .

In Figure 2.3, we show the fluid density $\rho$ , pressure $p$ , Mach number ${\left|\bm{u}\right|}/c_{s}$ (where $c_{s}=\sqrt{\gamma p/\rho}$ is the speed of sound), and magnetic pressure ${\left|\bm{b}\right|}^{2}/2$ computed on a uniform $200\times 200$ mesh at time $t=0.295$ . We note that our results are in good agreement with those reported in, e.g., [56, 42, 40]. In addition, it is emphasized in [5, 42, 56] that, due to rapid changes at the center of the rotation, many numerical methods produce oscillations or negative pressure values. We stress that during numerical simulations, we have not observed any oscillations, and the proposed PCCU scheme has produced no negative values of the computed pressure. The oscillation-free feature is further illustrated in Figure 2.4, where we zoom in on the center of the Mach number plots produced on three consecutively refined grids.

Example 4—Blast Problem.

In this example, we consider the blast problem first introduced in [5]. This benchmark problem is seen in a number of studies (see, e.g., [40, 39, 42, 65]) and is considered a challenge due to the low gas pressure and strong magnetosonic shocks. Negative pressures are easily produced near the shocks; see [39, 40] and references therein. The initial conditions are

(\rho,u,v,w,b_{1},b_{2},b_{3})(x,y,0)=\left(1,0,0,0,\frac{100}{\sqrt{4\pi}},0,0\right),\quad p(x,y,0)=\begin{cases}1000,&\sqrt{x^{2}+y^{2}}<0.1,\\ 0.1,&\mbox{othrwise}.\end{cases}

We take the specific heat ratio $\gamma=1.4$ and use zero-order extrapolation on the boundaries of the computational domain $[-0.5,0.5]\times[-0.5,0.5]$ .

The fluid density $\rho$ , pressure $p$ , magnitude of velocity ${\left|\bm{u}\right|}$ , and magnetic pressure ${\left|\bm{b}\right|}^{2}/2$ computed by the proposed PCCU scheme on a $200\times 200$ uniform mesh at $t=0.01$ are depicted in Figure 2.5. Additionally, the numerical experimentation of the proposed method resulted in positive pressure values throughout the entire computational domain, returning a minimum pressure of 0.10. Positive pressure values are also completely maintained when running the blast problem on a refined $400\times 400$ uniform grid (the fine mesh results are not shown here for brevity).

3 Shallow Water MHD

3.1 Governing Equations

In this section, we study the 2-D modified Godunov-Powell shallow water MHD system, which reads as

		$\displaystyle h_{t}+(hu)_{x}+(hv)_{y}=0,$		(3.1)
		$\displaystyle(hu)_{t}+\left(hu^{2}+\frac{g}{2}h^{2}-ha^{2}\right)_{x}+(huv-hab)_{y}=-a\left[(ha)_{x}+(hb)_{y}\right],$
		$\displaystyle(hv)_{t}+(huv-hab)_{x}+\left(hv^{2}+\frac{g}{2}h^{2}-hb^{2}\right)_{y}=-b\left[(ha)_{x}+(hb)_{y}\right],$
		$\displaystyle(ha)_{t}+(hbu-hav)_{y}=-u\left[(ha)_{x}+(hb)_{y}\right],$
		$\displaystyle(hb)_{t}+(hav-hbu)_{x}=-v\left[(ha)_{x}+(hb)_{y}\right].$

Here, $h$ is the fluid thickness, $u$ and $v$ represent the $x$ - and $y$ -velocity, $(a,b)^{\top}$ is the reduced magnetic field, which has units of velocity, and $g$ is the acceleration due to gravity. As in the ideal MHD system considered in §2, one can easily show that

(ha)_{x}+(hb)_{y}=0

(3.2)

as long as the field $(ha,hb)^{\top}$ is initially divergence-free. Therefore, the Godunov-Powell source terms on the RHS of (3.1) are theoretically zero. Still, they are added to the original shallow water MHD system (whose RHS is identically zero in the case of a flat bottom topography) to help enforce the divergence-free constraint (3.2) numerically; see, e.g., [23, 31, 58].

In order to develop a locally divergence-free numerical method for the system (3.1), this divergence constraint (3.2) must be enforced on the discrete level. As in §2, we achieve this goal by introducing the new variables $A:=(ha)_{x}$ and $B:=(hb)_{y}$ , differentiating the induction equation in (3.1) with respect to $x$ and $y$ , and obtaining the following equations for $A$ and $B$ ,

		$\displaystyle A_{t}+\big{(}uA-hbu_{y}\big{)}_{x}+\big{(}vA+hav_{x}\big{)}_{y}=0,$		(3.3)
		$\displaystyle B_{t}+\big{(}uB+hbu_{y}\big{)}_{x}+\big{(}vB-hav_{x}\big{)}_{y}=0.$		(3.3)

which are then added to the studied system (3.1).

Prior to introducing the numerical method for the augmented system (3.1), (3.3), we write it in the vector form

\bm{U}_{t}+\bm{F}(\bm{U})_{x}+\bm{G}(\bm{U})_{y}=Q(\bm{U})\bm{U}_{x}+R(\bm{U})\bm{U}_{y},

(3.4)

where $\bm{U}:=(h,hu,hv,ha,hb,A,B)^{\top}$ ,

	$\displaystyle\bm{F}(\bm{U})=\begin{pmatrix}hu\\ hu^{2}+\frac{g}{2}gh^{2}-ha^{2}\\ huv-hab\\ 0\\ hbu-hav\\ uA-hbu_{y}\\ uB+hbu_{y}\end{pmatrix},\quad Q(\bm{U})=\begin{pmatrix}0&0&0&0&0&0&0\\ 0&0&0&-a&0&0&0\\ 0&0&0&-b&0&0&0\\ 0&0&0&-u&0&0&0\\ 0&0&0&-v&0&0&0\\ 0&0&0&0&0&0&0\\ 0&0&0&0&0&0&0\end{pmatrix},$
	$\displaystyle\bm{G}(\bm{U})=\begin{pmatrix}hv\\ huv-hab\\ hv^{2}+\frac{g}{2}h^{2}-hb^{2}\\ hav-hbu\\ 0\\ vA+hav_{x}\\ vB-hav_{x},\end{pmatrix},\quad R(\bm{U})=\begin{pmatrix}0&0&0&0&0&0&0\\ 0&0&0&0&-a&0&0\\ 0&0&0&0&-b&0&0\\ 0&0&0&0&-u&0&0\\ 0&0&0&0&-v&0&0\\ 0&0&0&0&0&0&0\\ 0&0&0&0&0&0&0\end{pmatrix}.$

3.2 Numerical Method

We now extend the PCCU scheme developed in §2.2 to the shallow water MHD system.

Following the notation from section 2.2, the semi-discrete PCCU scheme still reads as (2.8)–(2.9), and the resulting system of ODEs is to be numerically integrated using an appropriate ODE solver, for instance, the three-stage third-order SSP Runge-Kutta, which we have used in the numerical experiments reported in §3.3.

In §3.2.1, §3.2.2, and §3.2.3 below, we focus on details of the scheme, which are different from the ideal MHD case.

3.2.1 Piecewise Linear Reconstruction

A piecewise linear reconstruction is performed for the discrete variables

\bm{W}_{j,k}:=(\,\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$h$\kern 0.0pt}}}_{j,k},u_{j,k},v_{j,k},(\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$ha$\kern 0.0pt}}})_{j,k},(\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$hb$\kern 0.0pt}}})_{j,k},\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$A$\kern 0.0pt}}}_{j,k},\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$B$\kern 0.0pt}}}_{j,k})^{\top},

where $u_{j,k}:=\nicefrac{{(\,\hbox{\vbox{\hrule height=0.5pt\kern 1.20552pt\hbox{\kern-0.35002pt$hu$\kern 0.0pt}}})_{j,k}}}{{\,\hbox{\vbox{\hrule height=0.5pt\kern 1.20552pt\hbox{\kern-0.35002pt$h$\kern 0.0pt}}}_{j,k}}}$ and $v_{j,k}:=\nicefrac{{(\,\hbox{\vbox{\hrule height=0.5pt\kern 1.20552pt\hbox{\kern-0.35002pt$hv$\kern 0.0pt}}})_{j,k}}}{{\,\hbox{\vbox{\hrule height=0.5pt\kern 1.20552pt\hbox{\kern-0.35002pt$h$\kern 0.0pt}}}_{j,k}}}$ . We then calculate cell interface values $\bm{W}^{\rm E,\rm W,\rm N,\rm S}_{j,k}$ using (2.11). The slopes $(W^{(i)}_{x})_{j,k}$ for $i\neq 4$ and $(W^{(i)}_{y})_{j,k}$ for $i\neq 5$ are computed using the generalized minmod limiter as in (2.12)–(2.13).

Like in the ideal MHD case, the slopes $((ha)_{x})_{j,k}$ and $((hb)_{y})_{j,k}$ are computed in a way that allows one to enforce the local discrete divergence-free condition $((ha)_{x})_{j,k}+((hb)_{y})_{j,k}\equiv 0$ for all $j,k$ . To this end, we proceed as in §2.2.1 and set

((ha)_{x})_{j,k}=\sigma_{j,k}\,\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$A$\kern 0.0pt}}}_{j,k},\quad((hb)_{y})_{j,k}=\sigma_{j,k}\,\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$B$\kern 0.0pt}}}_{j,k},

where

\sigma_{j,k}=\min\Big{\{}1,\sigma^{x}_{j,k},\sigma^{y}_{j,k}\Big{\}},

		$\displaystyle\sigma^{x}_{j,k}:=\left\{\begin{aligned} &\min\bigg{\{}1,\frac{((\widehat{ha})_{x})_{j,k}}{\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$A$\kern 0.0pt}}}_{j,k}}\bigg{\}}&&\mbox{if}\leavevmode\nobreak\ ((\widehat{ha})_{x})_{j,k}\,\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$A$\kern 0.0pt}}}_{j,k}>0,\\ &0&&\mbox{otherwise},\end{aligned}\right.$
		$\displaystyle\sigma^{y}_{j,k}:=\left\{\begin{aligned} &\min\bigg{\{}1,\frac{((\widehat{hb})_{y})_{j,k}}{\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$B$\kern 0.0pt}}}_{j,k}}\bigg{\}}&&\mbox{if}\leavevmode\nobreak\ ((\widehat{hb})_{y})_{j,k}\,\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$B$\kern 0.0pt}}}_{j,k}>0,\\ &0&&\mbox{otherwise},\end{aligned}\right.$

and $((\widehat{ha})_{x})_{j,k}$ and $((\widehat{hb})_{y})_{j,k}$ are computed using the generalized minmod reconstruction as in (2.12)–(2.13).

3.2.2 One-Sided Speeds of Propagation

We point out that in the shallow water MHD case, computing the one-sided speeds $s^{\pm}_{{j+\frac{1}{2}},k}$ and $s^{\pm}_{j,{k+\frac{1}{2}}}$ needed in the semi-discretization (2.8)–(2.9), is significantly easier than in the ideal MHD case. We follow the general recipe and estimate the $x$ - and $y$ -directional speeds using the largest and smallest eigenvalues of the matrices $\frac{\partial\bm{F}}{\partial\bm{U}}(\bm{U})-Q(\bm{U})$ and $\frac{\partial\bm{G}}{\partial\bm{U}}(\bm{U})-R(\bm{U})$ , respectively. This results in

		$\displaystyle s_{{j+\frac{1}{2}},k}^{+}=\max\left\{u_{j,k}^{\rm E}+\sqrt{\big{(}a_{j,k}^{\rm E}\big{)}^{2}+gh_{j,k}^{\rm E}},\,u_{j+1,k}^{\rm W}+\sqrt{\big{(}a_{j+1,k}^{\rm W}\big{)}^{2}+gh_{j+1,k}^{\rm W}},\,0\right\},$
		$\displaystyle s_{{j+\frac{1}{2}},k}^{-}=\min\left\{u_{j,k}^{\rm E}-\sqrt{\big{(}a_{j,k}^{\rm E}\big{)}^{2}+gh_{j,k}^{\rm E}},\,u_{j+1,k}^{\rm W}-\sqrt{\big{(}a_{j+1,k}^{\rm W}\big{)}^{2}+gh_{j+1,k}^{\rm W}},\,0\right\},$
		$\displaystyle s_{j,{k+\frac{1}{2}}}^{+}=\max\left\{v_{j,k}^{\rm N}+\sqrt{\big{(}b_{j,k}^{\rm N}\big{)}^{2}+gh_{j,k}^{\rm N}},\,v_{j,k+1}^{\rm S}+\sqrt{\big{(}b_{j,k+1}^{\rm S}\big{)}^{2}+gh_{j,k+1}^{\rm S}},\,0\right\},$
		$\displaystyle s_{j,{k+\frac{1}{2}}}^{-}=\min\left\{v_{j,k}^{\rm N}-\sqrt{\big{(}b_{j,k}^{\rm N}\big{)}^{2}+gh_{j,k}^{\rm N}},\,v_{j,k+1}^{\rm S}-\sqrt{\big{(}b_{j,k+1}^{\rm S}\big{)}^{2}+gh_{j,k+1}^{\rm S}},\,0\right\}.$

3.2.3 Discretization of Nonconservative Products

In order to evaluate the contribution of the nonconservative terms $Q(\bm{U})\bm{U}_{x}$ appearing on the RHS of (3.4), we again follow the lines of [11, 8] and evaluate the corresponding integrals exactly:

	$\displaystyle Q^{(1)}_{j,k}=Q^{(6)}_{j,k}=Q^{(7)}_{j,k}=Q^{(1)}_{\bm{\Psi},{j+\frac{1}{2}},k}=Q^{(6)}_{\bm{\Psi},{j+\frac{1}{2}},k}=Q^{(7)}_{\bm{\Psi},{j+\frac{1}{2}},k}=0,$
	$\displaystyle Q^{(2)}_{j,k}=-\int\limits_{x_{j-\frac{1}{2}}}^{x_{j+\frac{1}{2}}}\frac{\widetilde{ha}(x,y_{k})}{\widetilde{h}(x,y_{k})}((ha)_{x})_{j,k}\,{\rm d}x$
	$\displaystyle\hskip 19.91684pt=\left\{\begin{aligned} &-a_{j,k}\sigma_{j,k}\,\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$A$\kern 0.0pt}}}_{j,k}\Delta x&&\mbox{if}\leavevmode\nobreak\ (h_{x})_{j,k}=0,\\ &-\sigma_{j,k}\,\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$A$\kern 0.0pt}}}_{j,k}\left(\frac{(\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$ha$\kern 0.0pt}}})_{j,k}(h_{x})_{j,k}-\,\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$h$\kern 0.0pt}}}_{j,k}\sigma_{j,k}\,\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$A$\kern 0.0pt}}}_{j,k}}{((h_{x})_{j,k})^{2}}\,\ln\bigg{(}\frac{h^{\rm E}_{j,k}}{h^{\rm W}_{j,k}}\bigg{)}+\frac{\sigma_{j,k}\,\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$A$\kern 0.0pt}}}_{j,k}\Delta x}{(h_{x})_{j,k}}\right)&&\mbox{otherwise},\end{aligned}\right.$
	$\displaystyle Q^{(3)}_{j,k}=-\int\limits_{x_{j-\frac{1}{2}}}^{x_{j+\frac{1}{2}}}\frac{\widetilde{hb}(x,y_{k})}{\widetilde{h}(x,y_{k})}((ha)_{x})_{j,k}\,{\rm d}x$
	$\displaystyle\hskip 19.91684pt=\left\{\begin{aligned} &-b_{j,k}\sigma_{j,k}\,\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$A$\kern 0.0pt}}}_{j,k}\Delta x&&\mbox{if}\leavevmode\nobreak\ (h_{x})_{j,k}=0,\\ &-\sigma_{j,k}\,\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$A$\kern 0.0pt}}}_{j,k}\left(\frac{(\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$hb$\kern 0.0pt}}})_{j,k}(h_{x})_{j,k}-\,\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$h$\kern 0.0pt}}}_{j,k}((hb)_{x})_{j,k}}{((h_{x})_{j,k})^{2}}\,\ln\bigg{(}\frac{h^{\rm E}_{j,k}}{h^{\rm W}_{j,k}}\bigg{)}+\frac{((hb)_{x})_{j,k}\Delta x}{(h_{x})_{j,k}}\right)&&\mbox{otherwise},\end{aligned}\right.$
	$\displaystyle Q^{(4)}_{j,k}=-\int\limits_{x_{j-\frac{1}{2}}}^{x_{j+\frac{1}{2}}}\widetilde{u}(x,y_{k})((ha)_{x})_{j,k}\,{\rm d}x=-u_{j,k}\sigma_{j,k}\,\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$A$\kern 0.0pt}}}_{j,k}\Delta x,$
	$\displaystyle Q^{(5)}_{j,k}=-\int\limits_{x_{j-\frac{1}{2}}}^{x_{j+\frac{1}{2}}}\widetilde{v}(x,y_{k})((ha)_{x})_{j,k}\,{\rm d}x=-v_{j,k}\sigma_{j,k}\,\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$A$\kern 0.0pt}}}_{j,k}\Delta x,$
	$\displaystyle Q^{(2)}_{\bm{\Psi},{j+\frac{1}{2}},k}=-\int\limits_{0}^{1}\frac{(ha)^{\rm E}_{j,k}+s{\left((ha)^{\rm W}_{j+1,k}-(ha)^{\rm E}_{j,k}\right)}}{h^{\rm E}_{j,k}+s{\left(h^{\rm W}_{j+1,k}-h^{\rm E}_{j,k}\right)}}[ha]_{{j+\frac{1}{2}},k}\,{\rm d}s$
	$\displaystyle\hskip 19.91684pt=\left\{\begin{aligned} &-\frac{1}{2}\left(a^{\rm E}_{j,k}+a^{\rm W}_{j+1,k}\right)[ha]_{{j+\frac{1}{2}},k}&&\mbox{if}\leavevmode\nobreak\ [h]_{{j+\frac{1}{2}},k}=0,\\ &-[ha]_{{j+\frac{1}{2}},k}\left(\frac{(ha)^{\rm E}_{j,k}[h]_{{j+\frac{1}{2}},k}-h^{\rm E}_{j,k}[ha]_{{j+\frac{1}{2}},k}}{[h]_{{j+\frac{1}{2}},k}^{2}}\,\ln\bigg{(}\frac{h^{\rm W}_{j+1,k}}{h^{\rm E}_{j,k}}\bigg{)}+\frac{[ha]_{{j+\frac{1}{2}},k}}{[h]_{{j+\frac{1}{2}},k}}\right)&&\mbox{otherwise},\end{aligned}\right.$
	$\displaystyle Q^{(3)}_{\bm{\Psi},{j+\frac{1}{2}},k}=-\int\limits_{0}^{1}\frac{(hb)^{\rm E}_{j,k}+s{\left((hb)^{\rm W}_{j+1,k}-(hb)^{\rm E}_{j,k}\right)}}{h^{\rm E}_{j,k}+s{\left(h^{\rm W}_{j+1,k}-h^{\rm E}_{j,k}\right)}}[ha]_{{j+\frac{1}{2}},k}\,{\rm d}s$
	$\displaystyle\hskip 19.91684pt=\left\{\begin{aligned} &-\frac{1}{2}\left(b^{\rm E}_{j,k}+b^{\rm W}_{j+1,k}\right)[ha]_{{j+\frac{1}{2}},k}&&\mbox{if}\leavevmode\nobreak\ [h]_{{j+\frac{1}{2}},k}=0,\\ &-[ha]_{{j+\frac{1}{2}},k}\left(\frac{(hb)^{\rm E}_{j,k}[h]_{{j+\frac{1}{2}},k}-h^{\rm E}_{j,k}[hb]_{{j+\frac{1}{2}},k}}{[h]_{{j+\frac{1}{2}},k}^{2}}\,\ln\bigg{(}\frac{h^{\rm W}_{j+1,k}}{h^{\rm E}_{j,k}}\bigg{)}+\frac{[hb]_{{j+\frac{1}{2}},k}}{[h]_{{j+\frac{1}{2}},k}}\right)&&\mbox{otherwise},\end{aligned}\right.$
	$\displaystyle Q^{(4)}_{\bm{\Psi},{j+\frac{1}{2}},k}=-\int\limits_{0}^{1}\left\{u^{\rm E}_{j,k}+s{\left(u^{\rm W}_{j+1,k}-u^{\rm E}_{j,k}\right)}\right\}[ha]_{{j+\frac{1}{2}},k}\,{\rm d}s=-\frac{1}{2}\left(u^{\rm E}_{j,k}+u^{\rm W}_{j+1,k}\right)[ha]_{{j+\frac{1}{2}},k},$
	$\displaystyle Q^{(5)}_{\bm{\Psi},{j+\frac{1}{2}},k}=-\int\limits_{0}^{1}\left\{v^{\rm E}_{j,k}+s{\left(v^{\rm W}_{j+1,k}-v^{\rm E}_{j,k}\right)}\right\}[ha]_{{j+\frac{1}{2}},k}\,{\rm d}s=-\frac{1}{2}\left(v^{\rm E}_{j,k}+v^{\rm W}_{j+1,k}\right)[ha]_{{j+\frac{1}{2}},k},$

where

		$\displaystyle a_{j,k}:=\frac{(\,\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$ha$\kern 0.0pt}}})_{j,k}}{\,\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$h$\kern 0.0pt}}}_{j,k}},\quad b_{j,k}:=\frac{(\,\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$hb$\kern 0.0pt}}})_{j,k}}{\,\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$h$\kern 0.0pt}}}_{j,k}},\quad a^{\rm E(W)}_{j,k}:=\frac{(ha)^{\rm E(W)}_{j,k}}{h^{\rm E(W)}_{j,k}},\quad b^{\rm E(W)}_{j,k}:=\frac{(hb)^{\rm E(W)}_{j,k}}{h^{\rm E(W)}_{j,k}},$
		$\displaystyle[h]_{{j+\frac{1}{2}},k}:=h^{\rm W}_{j+1,k}-h^{\rm E}_{j,k},\quad[ha]_{{j+\frac{1}{2}},k}:=(ha)^{\rm W}_{j+1,k}-(ha)^{\rm E}_{j,k},\quad[hb]_{{j+\frac{1}{2}},k}:=(hb)^{\rm W}_{j+1,k}-(hb)^{\rm E}_{j,k}.$

The contribution of the nonconservative terms $R(\bm{U})\bm{U}_{x}$ appearing on the RHS of (3.4) is obtained in a similar manner and given by

	$\displaystyle R^{(1)}_{j,k}=R^{(6)}_{j,k}=R^{(7)}_{j,k}=R^{(1)}_{\bm{\Psi},j,{k+\frac{1}{2}}}=R^{(6)}_{\bm{\Psi},j,{k+\frac{1}{2}}}=R^{(7)}_{\bm{\Psi},j,{k+\frac{1}{2}}}=0,$
	$\displaystyle R^{(2)}_{j,k}=\left\{\begin{aligned} &-a_{j,k}\sigma_{j,k}\,\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$B$\kern 0.0pt}}}_{j,k}\Delta y&&\mbox{if}\leavevmode\nobreak\ (h_{y})_{j,k}=0,\\ &-\sigma_{j,k}\,\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$B$\kern 0.0pt}}}_{j,k}\left(\frac{(\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$ha$\kern 0.0pt}}})_{j,k}(h_{y})_{j,k}-\,\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$h$\kern 0.0pt}}}_{j,k}((ha)_{y})_{j,k}}{((h_{y})_{j,k})^{2}}\,\ln\bigg{(}\frac{h^{\rm N}_{j,k}}{h^{\rm S}_{j,k}}\bigg{)}+\frac{((ha)_{y})_{j,k}\Delta y}{(h_{y})_{j,k}}\right)&&\mbox{otherwise},\end{aligned}\right.$
	$\displaystyle R^{(3)}_{j,k}=\left\{\begin{aligned} &-b_{j,k}\sigma_{j,k}\,\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$B$\kern 0.0pt}}}_{j,k}\Delta y&&\mbox{if}\leavevmode\nobreak\ (h_{y})_{j,k}=0,\\ &-\sigma_{j,k}\,\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$B$\kern 0.0pt}}}_{j,k}\left(\frac{(\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$hb$\kern 0.0pt}}})_{j,k}(h_{y})_{j,k}-\,\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$h$\kern 0.0pt}}}_{j,k}\sigma_{j,k}\,\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$B$\kern 0.0pt}}}_{j,k}}{((h_{y})_{j,k})^{2}}\,\ln\bigg{(}\frac{h^{\rm N}_{j,k}}{h^{\rm S}_{j,k}}\bigg{)}+\frac{\sigma_{j,k}\,\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$B$\kern 0.0pt}}}_{j,k}\Delta y}{(h_{y})_{j,k}}\right)&&\mbox{otherwise},\end{aligned}\right.$
	$\displaystyle R^{(4)}_{j,k}=-u_{j,k}\sigma_{j,k}\,\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$B$\kern 0.0pt}}}_{j,k}\Delta y,\quad R^{(5)}_{j,k}=-v_{j,k}\sigma_{j,k}\,\hbox{\vbox{\hrule height=0.5pt\kern 1.72218pt\hbox{\kern-0.50003pt$B$\kern 0.0pt}}}_{j,k}\Delta y,$
	$\displaystyle R^{(2)}_{\bm{\Psi},j,{k+\frac{1}{2}}}=\left\{\begin{aligned} &-\frac{1}{2}\left(a^{\rm N}_{j,k}+a^{\rm S}_{j,k+1}\right)[hb]_{j,{k+\frac{1}{2}}}&&\mbox{if}\leavevmode\nobreak\ [h]_{j,{k+\frac{1}{2}}}=0,\\ &-[hb]_{j,{k+\frac{1}{2}}}\left(\frac{(ha)^{\rm N}_{j,k}[h]_{j,{k+\frac{1}{2}}}-h^{\rm N}_{j,k}[ha]_{j,{k+\frac{1}{2}}}}{[h]_{j,{k+\frac{1}{2}}}^{2}}\,\ln\bigg{(}\frac{h^{\rm S}_{j,k+1}}{h^{\rm N}_{j,k}}\bigg{)}+\frac{[ha]_{j,{k+\frac{1}{2}}}}{[h]_{j,{k+\frac{1}{2}}}}\right)&&\mbox{otherwise},\end{aligned}\right.$
	$\displaystyle R^{(3)}_{\bm{\Psi},j,{k+\frac{1}{2}}}=\left\{\begin{aligned} &-\frac{1}{2}\left(b^{\rm N}_{j,k}+b^{\rm S}_{j,k+1}\right)[hb]_{j,{k+\frac{1}{2}}}&&\mbox{if}\leavevmode\nobreak\ [h]_{j,{k+\frac{1}{2}}}=0,\\ &-[hb]_{j,{k+\frac{1}{2}}}\left(\frac{(hb)^{\rm N}_{j,k}[h]_{j,{k+\frac{1}{2}}}-h^{\rm N}_{j,k}[hb]_{j,{k+\frac{1}{2}}}}{[h]_{j,{k+\frac{1}{2}}}^{2}}\,\ln\bigg{(}\frac{h^{\rm S}_{j,k+1}}{h^{\rm N}_{j,k}}\bigg{)}+\frac{[hb]_{j,{k+\frac{1}{2}}}}{[h]_{j,{k+\frac{1}{2}}}}\right)&&\mbox{otherwise},\end{aligned}\right.$
	$\displaystyle R^{(4)}_{\bm{\Psi},j,{k+\frac{1}{2}}}=-\frac{1}{2}\left(u^{\rm N}_{j,k}+u^{\rm S}_{j,k+1}\right)[hb]_{j,{k+\frac{1}{2}}},\quad R^{(5)}_{\bm{\Psi},j,{k+\frac{1}{2}}}=-\frac{1}{2}\left(v^{\rm N}_{j,k}+v^{\rm S}_{j,k+1}\right)[hb]_{j,{k+\frac{1}{2}}}.$

where

		$\displaystyle a^{\rm N(S)}_{j,k}:=\frac{(ha)^{\rm N(S)}_{j,k}}{h^{\rm N(S)}_{j,k}},\quad b^{\rm N(S)}_{j,k}:=\frac{(hb)^{\rm N(S)}_{j,k}}{h^{\rm N(S)}_{j,k}},\quad[h]_{j,{k+\frac{1}{2}}}:=h^{\rm S}_{j,k+1}-h^{\rm N}_{j,k},$
		$\displaystyle[ha]_{j,{k+\frac{1}{2}}}:=(ha)^{\rm S}_{j,k+1}-(ha)^{\rm N}_{j,k},\quad[hb]_{j,{k+\frac{1}{2}}}:=(hb)^{\rm S}_{j,k+1}-(hb)^{\rm N}_{j,k}.$

3.3 Numerical Examples

In this section, we apply the proposed PCCU scheme to the 2-D shallow water MHD equations. In all of the examples, the CFL number is set to 0.25 and the minmod parameter is set to $\theta=1.3$ .

Example 5—Orszag-Tang-Like Problem.

This example taken from [18, 68] is similar to that of the ideal MHD Orszag-Tang problem studied in Example 2.

The shallow water MHD system is considered in the computational domain $[0,2\pi]\times[0,2\pi]$ subject to the periodic boundary conditions in both $x$ and $y$ -directions and the following smooth initial data:

(h,u,v,a,b)(x,y,0)=\left(\frac{25}{9},-\sin y,\sin x,-\sin y,\sin(2x)\right).

We compute the numerical solution by the proposed PCCU scheme on a uniform $200\times 200$ mesh until the final time $t=2$ . Time snapshots of $h$ and $\sqrt{a^{2}+b^{2}}$ at $t=1$ and 2 are plotted in Figure 3.1. As one can see, the initially smooth solution breaks down and develops multiple shock waves, whose interaction leads to the appearance of many essential features of MHD turbulence. We observe that the obtained results are in good agreement with those reported in [18, 68].

Example 6—Rotor-Like Problem.

Next, we consider a rotor-like problem taken from [18, 32]. This benchmark, which is an extension of the ideal MHD rotor problem studied in Example 3, portrays a disk with radius 0.1 of significant fluid depth $h$ rotating in a magnetic field.

The initial data

(h,u,v,ha,hb)=\begin{cases}(10,-y,x,1,0),&\sqrt{x^{2}+y^{2}}<0.1,\\ (1,0,0,1,0),&\mbox{otherwise},\end{cases}

are prescribed in the computational domain $[-1,1]\times[-1,1]$ and zero-order extrapolation boundary conditions are set along its boundary. The solution computed by the proposed PCCU scheme on a uniform $200\times 200$ mesh at time $t=0.2$ is plotted Figure 3.2. The obtained results are oscillation-free and overall comparable to those reported in [18, 32].

Example 7—Explosion Problem.

In the final example, we numerically solve the explosion problem studied in [32, 57]. This is another benchmark for the shallow water MHD equations considered subject to the following initial data:

(h,u,v,a,b)(x,y,0)=\begin{cases}(1,0,0,0.1,0),&\sqrt{x^{2}+y^{2}}<0.3,\\ (0.1,0,0,1,0),&\mbox{otherwise}.\end{cases}

In this example, we take the computational domain $[-1,1]\times[-1,1]$ and implement the zero-order extrapolation boundary conditions along its boundary.

The solution of the explosion problem consists of a shock traveling away from the center, a rarefaction wave traveling toward the origin, and two Alfvén waves. We compute the solution by the proposed PCCU scheme on a uniform $200\times 200$ mesh. The obtained results, shown in Figure 3.3 at $t=0.25$ , are non-oscillatory and agree well with the corresponding results in [32, 57].

4 Conclusion

In this paper, we have developed a new second-order unstaggered path-conservative central-upwind (PCCU) scheme for the ideal and shallow water magnetohydrodynamics (MHD) systems. The proposed scheme is (i) locally divergence-free; (ii) Riemann-problem-solver-free; (iii) high-resolution; (iv) robust; and (v) non-oscillatory. The derivation of the scheme is based on the Godunov-Powell nonconservative modifications of the studied MHD systems. The local divergence-free property is enforced by augmenting the studied systems with the evolution equations for the corresponding derivatives of the magnetic field components and by using these evolved quantities in the design of a special piecewise linear reconstruction of the magnetic field, which also guarantees a non-oscillatory nature of the resulting scheme. In addition, the proposed PCCU scheme allows for a proper treatment of the nonconservative product terms, which takes into account jumps of the normal component of the magnetic field across cell interfaces, thus providing stability. The performance of the new scheme has been illustrated on several benchmarks for both ideal and shallow water MHD systems producing high-resolution and oscillation-free results with positive computed quantities such as density, pressure, and water depth.

In our future work, we plan to develop a provably positivity-preserving high-order PCCU scheme as well as to introduce a new well-balanced PCCU scheme for more general shallow water MHD systems with the nonflat bottom topography and Coriolis forces are taken into account.

Acknowledgments

The work of A. Chertock and M. Redle were partially supported by NSF grants DMS-1818684 and DMS-2208438. The work of A. Kurganov was partially supported by NSFC grants 12111530004 and 12171226, and by the fund of the Guangdong Provincial Key Laboratory of Computational Science and Material Design (No. 2019B030301001). The work of K. Wu was partially supported by the NSFC grant 12171227. The authors would like to express their gratitude to Prof. Vladimir Zeitlin from the Laboratory of Dynamical Meteorology, Sorbonne University, Ecole Normale Supérieure, CNRS, Paris, for motivating and fruitful discussions.

References

[1] S. Ahmed and S. Zia, The higher-order CESE method for two-dimensional shallow water magnetohydrodynamics equations, Eur. J. Pure Appl. Math., 12 (2019), pp. 1464–1482.
[2] D. S. Balsara, Divergence-free reconstruction of magnetic fields and WENO schemes for magnetohydrodynamics, J. Comput. Phys., 228 (2009), pp. 5040–5056.
[3] D. S. Balsara, R. Kumar, and P. Chandrashekar, Globally divergence-free DG scheme for ideal compressible MHD, Commun. Appl. Math. Comput. Sci., 16 (2021), pp. 59–98.
[4] D. S. Balsara, T. Rumpf, M. Dumbser, and C.-D. Munz, Efficient, high accuracy ADER-WENO schemes for hydrodynamics and divergence-free magnetohydrodynamics, J. Comput. Phys., 228 (2009), pp. 2480–2516.
[5] D. S. Balsara and D. S. Spicer, A staggered mesh algorithm using high order Godunov fluxes to ensure solenoidal magnetic fields in magnetohydrodynamic simulations, J. Comput. Phys., 149 (1999), pp. 270–292.
[6] J. U. Brackbill and D. C. Barnes, The effect of nonzero $\nabla\cdot B$ on the numerical solution of the magnetohydrodynamic equations, J. Comput. Phys., 35 (1980), pp. 426–430.
[7] M. Brio and C. C. Wu, An upwind differencing scheme for the equations of ideal magnetohydrodynamics, J. Comput. Phys., 75 (1988), pp. 400–422.
[8] M. J. Castro Díaz, A. Kurganov, and T. Morales de Luna, Path-conservative central-upwind schemes for nonconservative hyperbolic systems, ESAIM Math. Model. Numer. Anal., 53 (2019), pp. 959–985.
[9] P. Chandrashekar and C. Klingenberg, Entropy stable finite volume scheme for ideal compressible mhd on 2-d cartesian meshes, SIAM J. Numer. Anal., 54 (2016), pp. 1313–1340.
[10] A. J. Christlieb, J. A. Rossmanith, and Q. Tang, Finite difference weighted essentially non-oscillatory schemes with constrained transport for ideal magnetohydrodynamics, J. Comput. Phys., 268 (2014), pp. 302–325.
[11] S. Chu, A. Kurganov, and M. Na, Fifth-order A-WENO schemes based on the path-conservative central-upwind method, J. Comput. Phys., 469 (2022). Paper No. 111508, 22 pp.
[12] W. Dai and P. R. Woodward, A simple finite difference scheme for multidimensional magnetohydrodynamical equations, J. Comput. Phys., 142 (1998), pp. 331–369.
[13] A. Dedner, F. Kemm, D. Kröner, C.-D. Munz, T. Schnitzer, and M. Wesenberg, Hyperbolic divergence cleaning for the mhd equations, J. Comput. Phys., 175 (2002), pp. 645–673.
[14] P. J. Dellar, A note on magnetic monopoles and the one-dimensional MHD Riemann problem, J. Comput. Phys., 172 (2001), pp. 392–398.
[15] , Dispersive shallow water magnetohydrodynamics, Phys. Plasmas, 10 (2003), pp. 581–590.
[16] D. Derigs, G. J. Gassner, S. Walch, and A. R. Winters, Entropy stable finite volume approximations for ideal magnetohydrodynamics, Jahresber. Dtsch. Math.-Ver., 120 (2018), pp. 153–219.
[17] C. R. DeVore, Flux-corrected transport techniques for multidimensional compressible magnetohydrodynamics, J. Comput. Phys., 92 (1991), pp. 142–160.
[18] J. Duan and H. Tang, High-order accurate entropy stable finite difference schemes for the shallow water magnetohydrodynamics, J. Comput. Phys., 431 (2021). Paper No. 110136, 26 pp.
[19] M. Dumbser, D. S. Balsara, M. Tavelli, and F. Fambri, A divergence-free semi-implicit finite volume scheme for ideal, viscous, and resistive magnetohydrodynamics, Internat. J. Numer. Methods Fluids, 89 (2019), pp. 16–42.
[20] B. Einfeldt, C.-D. Munz, P. L. Roe, and B. Sjögreen, On Godunov-type methods near low densities, J. Comput. Phys., 92 (1991), pp. 273–295.
[21] C. R. Evans and J. F. Hawley, Simulation of magnetohydrodynamic flows: A constrained transport method, Astrophys. J., 332 (1988), pp. 659–677.
[22] P. Fu, F. Li, and Y. Xu, Globally divergence-free discontinuous Galerkin methods for ideal magnetohydrodynamic equations, J. Sci. Comput., 77 (2018), pp. 1621–1659.
[23] F. G. Fuchs, A. D. McMurry, S. Mishra, N. H. Risebro, and K. Waagan, Approximate Riemann solvers and robust high-order finite volume schemes for multi-dimensional ideal MHD equations, Commun. Comput. Phys., 9 (2011), pp. 324–362.
[24] T. A. Gardiner and J. M. Stone, An unsplit Godunov method for ideal MHD via constrained transport, J. Comput. Phys., 205 (2005), pp. 509–539.
[25] P. A. Gilman, Magnetohydrodynamic “shallow water” equations for the solar tachocline, Astrophys. J. Lett., 544 (2000), pp. L79–L82.
[26] S. K. Godunov, Symmetric form of the equations of magnetohydrodynamics (in Russian), Numerical Methods for Mechanics of Continuum Medium, 1 (1972), pp. 26–34.
[27] S. Gottlieb, D. Ketcheson, and C.-W. Shu, Strong stability preserving Runge-Kutta and multistep time discretizations, World Scientific Publishing Co. Pte. Ltd., Hackensack, NJ, 2011.
[28] S. Gottlieb, C.-W. Shu, and E. Tadmor, Strong stability-preserving high-order time discretization methods, SIAM Rev., 43 (2001), pp. 89–112.
[29] C. Helzel, J. A. Rossmanith, and B. Taetz, An unstaggered constrained transport method for the 3D ideal magnetohydrodynamic equations, J. Comput. Phys., 230 (2011), pp. 3803–3829.
[30] , A high-order unstaggered constrained-transport method for the three-dimensional ideal magnetohydrodynamic equations based on the method of lines, SIAM J. Sci. Comput., 35 (2013), pp. A623–A651.
[31] P. Janhunen, A positive conservative method for magnetohydrodynamics based on HLL and Roe methods, J. Comput. Phys., 160 (2000), pp. 649–661.
[32] T. Kröger and M. Lukáčová-Medviďová, An evolution Galerkin scheme for the shallow water magnetohydrodynamic equations in two space dimensions, J. Comput. Phys., 206 (2005), pp. 122–149.
[33] A. Kurganov and C.-T. Lin, On the reduction of numerical dissipation in central-upwind schemes, Commun. Comput. Phys., 2 (2007), pp. 141–163.
[34] A. Kurganov, S. Noelle, and G. Petrova, Semidiscrete central-upwind schemes for hyperbolic conservation laws and Hamilton-Jacobi equations, SIAM J. Sci. Comput., 23 (2001), pp. 707–740.
[35] A. Kurganov and E. Tadmor, New high-resolution central schemes for nonlinear conservation laws and convection-diffusion equations, J. Comput. Phys., 160 (2000), pp. 241–282.
[36] , Solution of two-dimensional Riemann problems for gas dynamics without Riemann problem solvers, Numer. Methods Partial Differential Equations, 18 (2002), pp. 584–608.
[37] N. Lahaye and V. Zeitlin, Coherent magnetic modon solutions in quasi-geostrophic shallow water magnetohydrodynamics, J. Fluid Mech., 941 (2022).
[38] F. Li and C.-W. Shu, Locally divergence-free discontinuous Galerkin methods for MHD equations, J. Sci. Comput., 22/23 (2005), pp. 413–442.
[39] F. Li and L. Xu, Arbitrary order exactly divergence-free central discontinuous Galerkin methods for ideal MHD equations, J. Comput. Phys., 231 (2012), pp. 2655–2675.
[40] F. Li, L. Xu, and S. Yakovlev, Central discontinuous Galerkin methods for ideal MHD equations with the exactly divergence-free magnetic field, J. Comput. Phys., 230 (2011), pp. 4828–4847.
[41] K.-A. Lie and S. Noelle, On the artificial compression method for second-order nonoscillatory central difference schemes for systems of conservation laws, SIAM J. Sci. Comput., 24 (2003), pp. 1157–1174.
[42] M. Liu, M. Zhang, C. Li, and F. Shen, A new locally divergence-free WLS-ENO scheme based on the positivity-preserving finite volume method for ideal MHD equations, J. Comput. Phys., 447 (2021). Paper No. 110694, 24 pp.
[43] Y. Liu, C.-W. Shu, and M. Zhang, Entropy stable high order discontinuous Galerkin methods for ideal compressible MHD on structured meshes, J. Comput. Phys., 354 (2018), pp. 163–178.
[44] P. Londrillo and L. Del Zanna, On the divergence-free condition in Godunov-type schemes for ideal magnetohydrodynamics: the upwind constrained transport method, J. Comput. Phys., 195 (2004), pp. 17–48.
[45] S. Mishra and E. Tadmor, Constraint preserving schemes using potential-based fluxes. III. Genuinely multi-dimensional schemes for MHD equations, ESAIM Math. Model. Numer. Anal., 46 (2012), pp. 661–680.
[46] H. Nessyahu and E. Tadmor, Nonoscillatory central differencing for hyperbolic conservation laws, J. Comput. Phys., 87 (1990), pp. 408–463.
[47] S. A. Orszag and C.-M. Tang, Small-scale structure of two-dimensional magnetohydrodynamic turbulence, J. Fluid Mech., 90 (1979), pp. 129–143.
[48] A. Petrosyan, D. Klimachkov, M. Fedotova, and T. Zinyakov, Shallow water magnetohydrodynamics in plasma astrophysics. Waves, turbulence, and zonal flows, Atmosphere, 11 (2020). Paper No. 314, 16 pp.
[49] K. G. Powell, An approximate Riemann solver for magnetohydrodynamics, in Upwind and High-Resolution Schemes, M. Y. Hussaini, B. van Leer, and J. Van Rosendale, eds., Springer Berlin Heidelberg, Berlin, Heidelberg, 1997, pp. 570–583.
[50] K. G. Powell, P. L. Roe, T. J. Linde, T. I. Gombosi, and D. L. De Zeeuw, A solution-adaptive upwind scheme for ideal magnetohydrodynamics, J. Comput. Phys., 154 (1999), pp. 284–309.
[51] K. G. Powell, P. L. Roe, R. S. Myong, T. Gombosi, and D. De Zeeuw, An upwind scheme for magnetohydrodynamics, in 12th Computational Fluid Dynamics Conference: AIAA Paper 95-1704-CP, 1995, pp. 661–674.
[52] S. Qamar and G. Warnecke, Application of space-time CE/SE method to shallow water magnetohydrodynamic equations, J. Comput. Appl. Math., 196 (2006), pp. 132–149.
[53] J. A. Rossmanith, An unstaggered, high-resolution constrained transport method for magnetohydrodynamic flows, SIAM J. Sci. Comput., 28 (2006), pp. 1766–1797.
[54] D. Ryu, F. Miniati, T. W. Jones, and A. Frank, A divergence-free upwind code for multidimensional magnetohydrodynamic flows, Astrophys. J., 509 (1998), pp. 244–255.
[55] P. K. Sweby, High resolution schemes using flux limiters for hyperbolic conservation laws, SIAM J. Numer. Anal., 21 (1984), pp. 995–1011.
[56] G. Tóth, The $\nabla\cdot B=0$ constraint in shock-capturing magnetohydrodynamics codes, J. Comput. Phys., 161 (2000), pp. 605–652.
[57] R. Touma, Unstaggered central schemes with constrained transport treatment for ideal and shallow water magnetohydrodynamics, Appl. Numer. Math., 60 (2010), pp. 752–766.
[58] K. Waagan, C. Federrath, and C. Klingenberg, A robust numerical scheme for highly compressible magnetohydrodynamics: nonlinear stability, implementation and tests, J. Comput. Phys., 230 (2011), pp. 3331–3351.
[59] A. R. Winters and G. J. Gassner, An entropy stable finite volume scheme for the equations of shallow water magnetohydrodynamics, J. Sci. Comput., 67 (2016), pp. 514–539.
[60] K. Wu, Positivity-preserving analysis of numerical schemes for ideal magnetohydrodynamics, SIAM J. Numer. Anal., 56 (2018), pp. 2124–2147.
[61] K. Wu and C.-W. Shu, Geometric quasilinearization framework for analysis and design of bound-preserving schemes, SIAM Review. To appear.
[62] , A provably positive discontinuous Galerkin method for multidimensional ideal magnetohydrodynamics, SIAM J. Sci. Comput., 40 (2018), pp. B1302–B1329.
[63] , Provably positive high-order schemes for ideal magnetohydrodynamics: analysis on general meshes, Numer. Math., 142 (2019), pp. 995–1047.
[64] Z. Xu, D. S. Balsara, and H. Du, Divergence-free WENO reconstruction-based finite volume scheme for solving ideal MHD equations on triangular meshes, Commun. Comput. Phys., 19 (2016), pp. 841–880.
[65] S. Yakovlev, L. Xu, and F. Li, Locally divergence-free central discontinuous Galerkin methods for ideal MHD equations, J. Comput. Sci., 4 (2013), pp. 80–91.
[66] V. Zeitlin, Remarks on rotating shallow-water magnetohydrodynamics, Nonlinear Process. Geophys., 20 (2013), pp. 893–898.
[67] V. Zeitlin, C. Lusso, and F. Bouchut, Geostrophic vs magneto-geostrophic adjustment and nonlinear magneto-inertia-gravity waves in rotating shallow water magnetohydrodynamics, Geophys. Astrophys. Fluid Dyn., 109 (2015), pp. 497–523.
[68] S. Zia, M. Ahmed, and S. Qamar, Numerical solution of shallow water magnetohydrodynamic equations with non-flat bottom topography, Int. J. Comput. Fluid Dyn., 28 (2014), pp. 56–75.