Skew-Symmetric Splitting and Stability of High Order Central Schemes

Skew-symmetric splittings of the inviscid flux derivative for high order central schemes are studied and developed to improve their numerical stability without added high order numerical dissipation for long time wave propagations and long time integration of compressible turbulent flows. For flows containing discontinuities and multiscale turbulence fluctuations the Yee & Sjogreen [33] and Kotov et al. [15, 14] high order nonlinear filter approach is utilized in conjunction with the skew-symmetric form of high order central schemes. Due to the incomplete hyperbolic nature of the conservative ideal magnetohydrodynamics (MHD) governing equations, not all of the skew-symmetric splittings for gas dynamics can be extended to the ideal MHD. For the MHD the Ducros et al. [6] variants are constructed. In addition, four formulations of the MHD are considered: (a) the conservative MHD, (b) the Godunov/Powell non-conservative form, (c) the Janhunen MHD with magnetic field source terms [13], and (d) a MHD with source terms of [3]. The different formulation of the equations in conjunction with the variants of Ducros et al. type skew-symmetric splitting will be shown to have a strong effect on the stability of non-dissipative approximations. Representative test cases for both smooth flows and problems containing discontinuities for the ideal MHD are included. The results illustrate the improved stability by using the skew-symmetric splitting as part of the central base scheme instead of the pure high order central scheme.


Introduction
Starting the early 80s skew-symmetric splitting of certain component of the inviscid flux derivatives in conjunction with central schemes was shown to help with numerical stability for long time integration. For a certain splitting it can provide a stable energy norm estimate for the Euler equations with smooth flows. For other skew-symmetric formulations they can provide a discrete momentum conservation or a discrete kinetic preservation property. See Arakawa, Blaisdell et [1,2,35,6,32,33,21,15,14] for some discussions and performance of the combined approach for DNS and LES applications. Some of the skew-symmetric splitting for the gas dynamics flux derivatives are not applicable and/or cannot straightforwardly be extended to the ideal MHD [35]. Their construction is also dependent on the MHD governing equation formulation.
Entropy conservative schemes [27,5,31] are another class of methods that might have better stability properties than straightforward pure centered discretizations. In view of the fact that methods proposed in [27,5,31] are low order and their linear numerical dissipation approaches for shock-capturing require further improvement, Sjögreen  to construct a form of the high order entropy conservative numerical fluxes. Starting with the high order entropy conservative development of Sjögreen & Yee [21] for gas dynamics smooth flows, construction of efficient high order conservative numerical fluxes for problems containing discontinuities and for the ideal MHD are reported in Sjögreen & Yee [23,24]. Again, extension of high order entropy conservative numerical fluxes that were developed for gas dynamics to the MHD is not straightforward due to the non-strictly hyperbolic nature of the conservative ideal MHD equations. See [23,24]. Some high order comparisons between skew-symmetric splittings with high order entropy conservative centered schemes can be found in [23,24,22].
Objective and Outline: The objective of this paper is to develop and test four variants of Ducros et al. [6] type skew-symmetric splitting methods for equations of compressible ideal MHD, including problems with discontinuities. Their stability for problems with smooth solution will be investigated. Putting the discrete div B error issue aside, strictly speaking, there are several slightly different, but equivalent, ways to formulate the equations of MHD. We will show that the exact formulation of the equations will have a strong effect on the stability of nondissipative approximations. Specifically, four formulations of the MHD are considered: (a) the conservative MHD, (b) the Godunov/Powell non-conservative form, (c) the Janhunen MHD with magnetic field source terms [13], and (d) the MHD source term of [3]. Another complication is that not all skew-symmetric splitting for the perfect gas dynamics can be extended to the MHD. Furthermore, extension of some of these split schemes for the MHD consist of many variants due to the incomplete hyperbolic nature of the conservative MHD equations. For the Ducros [15,14] high order nonlinear filter approach is extended to include the four skewsymmetric splittings for the MHD as part (pre-processing step) of the base scheme. Different splitting formulations and comparison among the aforementioned four forms of the MHD will be illustrated.
2. An overview of skew-symmetric split approximations for gas dynamics Standard centered difference approximations of nonlinear conservation laws normally encounter nonlinear instabilities after a short time integration without added numerical dissipation. It is well known that the appearance of these instabilities can be delayed if the convective flux derivatives are written in an equivalent desired split form before the pure central approximation is employed. Hereafter this is referred to as a split approximation.
For example, a split approximation starts from rewriting the derivative of the product (ab) x as before discretization. Here a and b are functions of x and α, γ and β are parameters so chosen to be still equivalent to the original (ab) x before discretization. A common split derivative is by setting α = γ = β = 1/2 resulting in the form These methods have a long history in finite difference approximations; see, .e.g., [1,16]. See also a generalized conservative split convective derivative operators study by Pirozzoli [19] key mathematical idea is that formulas of type (2) can be used to estimate the L 2 norm or the energy norm of the computed solution. From physical considerations some of the splittings provide the discrete conservation of momentum or preservation of discrete kinetic energy. A well-known example is the linear system of conservation laws where A(x) is a symmetric matrix, and we solve for the unknown vector u = u(x, t) from given initial data u(x, 0) = u 0 (x). Boundary data are given at x = 0 and x = L. To show how this is done, e.g., we rewrite (3) in an mathematically equivalent form: and define the scalar product and norm by A norm estimate is obtained if (4) is multiplied by u and integrated over [0, L]. We obtain where the second equality is obtained from partial integration of (u, (Au) x ), and from the symmetry of A which allows it to be moved between the arguments of the scalar product. If the boundary data are such that u T Au| L 0 ≥ 0, then the estimate 1 2 holds, which under the assumption that max x |A x | is bounded leads to a stability estimate by use of Gronwall's lemma. Let x j = j∆x, j = 0, . . . , N be a grid with spacing ∆x, and let u j (t) denote a numerical approximation of u(x j , t). Consider the semi-discrete approximation of (4) where D is a centered finite difference operator approximating d/dx. Note that A(x) is a given function, so that the exact derivative A x can be used in (8). The discrete scalar product and norm are defined by where ω j > 0 are weights that are equal to one at most grid points, but are given special values near the boundaries j = 0 and j = N . The boundary modified norm weights together with special boundary modifications of D lead to the summation-by-parts property, see [28] for details. Thanks to the summation-by-parts property, the same technique that led to the estimate (7) can be used to obtain the semi discrete estimate The possible growth rate is determined by A x in both (7) and (9), so that the discrete estimate will have the same growth rate as the estimate of the continuous problem.
Ducros et al. Type Conservative Splitting: For nonlinear systems, such as the Euler equations of gas dynamics, split approximations have been used for a long time see, e.g., Ducros et al and Blaisdel et al. [6,2]. The split approximations makes use of (2) to rewrite different terms in the Euler equations as sums of three terms. The terms of the split form (2) are approximated by where D is a finite difference operator, and a and b are functions of x. As shown in Ducros et al. [6], the approximation (10) can be written in conservation form. For example, with the second order operator Du j = (u j+1 − u j−1 )/(2∆x), it holds that where ∆ + q j = (q j+1 − q j ). Equation (11) can be generalized to standard centered difference operators of 2p-th order of accuracy, The coefficients α To derive the conservative form of the split approximation for an arbitrary operator, the right hand side of the algebraic identity is written on conservative form by The conservative form of the split approximation becomes where the numerical flux is defined by To simplify the formulas of the conservative form of split approximations for systems of equations, define For the three dimensional Euler equations of gas dynamics equations, the x-direction flux is where the inviscid flux f = f (u), the velocity vector u = (u(x, y, z), v(x, y, z), z(x, y, z)) T in the x-, y-, and z-directions, ρ denotes the density, p is the pressure, and e is the total energy. Denote u j = (u j , v j , w j ) T as the discretization at the j grid location with the y and z discretiztion indices suppress for simplicity. The flux components can be written as products of two factors in many different ways, leading to different split approximations. One Ducros et al. split-type approximation of the gas dynamics flux derivative that will be used in this study is given by which by (17) can be written on conservative form with numerical flux function The more compact notation introduced in (18) allows (20) to be rewritten as The homogeneity property of the perfect gas dynamics inviscid flux implies that f (u) = A(u)u, where A(u) is the Jacobian of f (u). To make use of the homogeneity property, a nonconservative natural splitting is where the discretization is Here A x and D p A denote element wise application of differentiation and differencing respectively. The approximation (23) can be rewritten of conservative form with numerical flux where A k,m denotes element (k, m) of the matrix valued function A(x), and u m denotes the mth componen of the vector u.
A Semi-conservative Entropy Splitting of the Euler Flux Derivatives: Another splitting that gives entropy stability of the Euler equation of gas dynamics is by Olsson & Oliger, Gerritsen & Olsson, and Yee et al. [17,11,35]. They made used of Harten's symmetrizable form of the Euler equations in terms of the entropy variables [9] to obtain a semidiscrete splitting of the Euler equations with a discrete entropy stability by the summation-byparts approach. During the computations, the entropy splitting is written in terms of the sum of a conservative portion for the interior scheme (interior grid points) and with a summationby-parts for the boundary scheme (boundary points). Note that the Harten [9] and Gerritsen & Olsson entropy splitting form selects the un-physical branch of the inequality and was later corrected by Yee et al., hereafter referred to as the entropy splitting of the Euler equations. It is considered to be a semi-conservative splitting except at the boundary grid points. The entropy splitting of Olsson & Oliger, Gerritsen & Olsson, and Yee et al. [17,11,35] is a splitting which is of a form that is more suitable for the discrete stable energy norm estimate technique, including boundary scheme estimate for arbitrary order of central spatial schemes. See Yee et al. [36] for the formulation. For the 1D Euler equations the inviscid flux derivative f (u) x for a perfect gas is split into the following via the entropy variables W discussed in Harten [9]. where and See Yee et al. [36,35,25] for the formulation, the choice for α and numerical examples.
Several split discretizations were compared in [10] where discretization by the entropy splitting form was shown by numerical experiments to be one of the best performing for smooth flows. For their skew-symmetric splitting extension to the ideal MHD, see Sjögreen

MHD formulations
The equations of magnetohydrodynamics (MHD) is the system of conservation laws where the unknown field vector is where ρ is density, (u, v, w) are the velocities in the x, y and z directions, e the internal energy and (B (x) , B (y) , B (z) ) the magnetic field components. The x-direction flux is given by and similar expressions hold for g and h. Here p denotes the pressure where w = (u, v, w) and B = (B (x) , B (y) , B (z) ). In the last term on the left hand side of (28) multiplies div B. This term could be removed, because div B = 0 from a physical standpoint. Numerically, if we can keep the discrete form of div B to be very small, e.g., in the order of 10 −6 , this non-conservative form is a symmetrizable system with a complete set of eigenvectors. Godunov [12] was the first to show that (28) can be symmetrized. Some variants of (28) are to replace e by either e J = (0, 0, 0, 0, 0, u, v, w) T as suggested in [13] or by T as suggested in [3]. The conservative form (28) where the fifth component is The flux derivative approximation (30) can, by (17), be written on conservative form with numerical flux function where and The more compact notation introduced in (18) allows (30) to be rewritten as The notation in (35) will be used to describe the remaining three split forms. The second MHD split scheme is similar, but the magnetic field terms are split without density weighting. This splitting is denoted by DS2, and its numerical flux is given by The third splitting does not split the magnetic field components. It is denoted by DS3, and have numerical flux The fourth splitting, denoted by DS4, differs from DS3 only in that the first term of the three The performance of (30)-(38) will be evaluated by numerical experiments in Section 6.

Skew-Symmetric Splitting as Preprocessing
Step in the Framework of the Nonlinear Filter Method of Yee & Sjögreen [33] To make the discussion more self-contained, this section gives a brief overview of the high order nonlinear filter scheme of Yee et al. and Yee & Sjögreen [35,34,32,33] for accurate computations of DNS and LES of compressible turbulence for a wide range of flow types by introducing as little shock-capturing numerical dissipation as possible. Preprocessing Step by Skew-symmetric Splitting for Gas Dynamics: Before the application of a high-order non-dissipative spatial base scheme, a preprocessing step is employed to improve numerical stability. The inviscid flux derivatives of the governing equations are split into the following two ways, depending on the flow types and the desire for rigorous mathematical analysis or physical argument.
• Entropy splitting of [35] or the natural splitting described Section 3: These are nonconservative splittings and they are among some of the best in improving numerical stability for non-dissipative central schemes, especially for long time integration of shockfree turbulence. • The Ducros et al. splitting [6] for systems (or variants of the conservative skew-symmetric splitting described earlier): These are conservative splitting and are suitable for problems with discontinuities.

Base Scheme
Step Using the Preprocessing Step: A full time step is advanced using a high order non-dissipative (or very low dissipation) spatially central scheme on the split form of the governing partial differential equations (PDEs) (i.e., after the preprocessing step).
For the base scheme step, a full time step of high-order temporal discretization such as the third-order or fourth-order Runge-Kutta (RK3 or RK4) method is used. It is remarked that other temporal discretizations can be used for the base scheme step.

Post-Processing (Nonlinear Filter
Step): To further improve the accuracy of the computed solution from the base scheme step, after a full time step of a non-dissipative high order spatial base scheme, the post-processing step is used to nonlinearly filter the solution by a dissipative portion of a high-order shock-capturing scheme with a local flow sensor. The flow sensor provides locations and amounts of built-in shock-capturing dissipation that can be further reduced or eliminated. At each grid point a local flow sensor is employed to analyze the regularity of the computed flow data. Only the strong discontinuity locations would receive the full amount of shock-capturing dissipation. In smooth regions no shock-capturing dissipation would be added,  [15,14] and Sjögreen & Yee [23,24] for more details and the formula for corresponding numerical fluxes for the filter step (filter numerical fluxes). Filter numerical fluxes are not to be confused with the standard numerical fluxes for the full shock-capturing method.
The filter numerical fluxes for the filter step only involve the inviscid flux derivatives regardless if the flow is viscous or inviscid. If viscous terms are present, a matching high order central difference operator (as the inviscid difference operator) is included on the base scheme step. In viscous flow with non-periodic boundaries, for ease of summation-by-parts (SBP) numerical boundary closure [18] implementation for the viscous flux derivatives, the same inviscid central difference operator for the first derivative is employed twice for the viscous flux derivatives. This is due to the fact that if higher than second-order central base schemes are used, proper numerical boundary conditions in the vicinity of the non-periodic boundaries are needed for numerical stability and accuracy consideration. SBP numerical boundary closures are designed to help numerical stability and accuracy of the overall high order central scheme.
Remark For the gas dynamics the post-processing (nonlinear filter step) are applied to all of the equation set. For the MHD in uniform Cartesian grid, in order to obtain zero discrete div B error without any div B cleaning, the nonlinear filter step does not employ for the three magnetic field equations. See Yee & Sjögreen [34] for details.
The aforementioned high order nonlinear filter method is valid for the four forms of the MHD formulation and the four skew-symmetric splitting of the MHD to be used as the preprocessing step. In addition, the aforementioned high order nonlinear filter method is valid for the four forms of the MHD formulation and the different high order entropy conservative numerical fluxes as the spatial base schemes discussed in Sections 4 and 5 of Sjögreen & Yee [23,24].

Numerical results
For smooth flow test cases, only the high order base schemes are activated. For problems containing shocks, the aforementioned nonlinear filter approach of Yee & Sjögreen is utilized. For all test cases we use the classical fourth-order Runge-Kutta time discretization. Here, for illustration purposes, only one smart flow sensor (among the many variants indicated in [33] and Kotov et al. [15,14]) is chosen for the numerical experiment for the nonlinear filter approach. It is the third-order B-spline wavelet flow sensor developed in Sjögreen & Yee [20]. Due to a page limitation, many details of the development and extensive numerical testing with more representative test cases among the different formulations and the different governing equation sets will be reported in [22] for journal publication. For 3D DNS and LES gas dynamics computations include shock-free turbulence, turbulence with moderate and strong shocks are included in [15,14,24].
We will consider a two dimensional wave, traveling in the direction n = (cos α, sin α, 0), where α = 30 • . The size of the computational domain is 0 ≤ x ≤ 1/ cos α, and 0 ≤ y ≤ 1/ sin α with periodic boundary conditions. For this computation, A = 0.1. The grid spacing used was ∆x = ∆y = 0.02. The magnetic field derivatives in the term e div B in (28), when present, are approximated by standard centered difference operators of the same order of accuracy as used to approximate the flux derivatives. Figure 1 shows the evolution of the norm of the error for the four split forms DS1, DS2, DS3 and DS4 (30)-(38) inconjunction with the fourth-order central scheme. The comparison also include the standard fourth-order pure centered unsplit approximation (C04). The error is computed as the maximum norm of the pointwise error over x, y, and the eight solution components. The unsplit scheme is plotted in blue, C04DS1 in black (C04 in conjunction with DS1 split form), C04DS2 in red (C04 in conjunction with DS2 split form), C04DS3 in green (C04 in conjunction with DS3 split form), and C04DS4 in cyan color (C04 inconjunction with DS4 split form). The top left subplot of Fig. 1 depicts results obtained when e = e 0 (the conservative equations), and the top right subplot shows e = e G , the bottom left subplot uses e = e J , and the right subplot uses e = e B . The unsplit scheme, C04, and the splitting with unsplit magnetic field variables, C04DS3, perform similarly for all three forms of the equations. The reason could be that the discretization of div B is zero up to round-off errors for these two schemes, so that the exact form of e has a negligible effect on the result. The other schemes C04DS1, C04DS2, and C04DS4 have some split terms in the magnetic field equations, which prevents perfect discrete div B conservation. Interestingly, the presence of e G (top right subplot  Fig. 1) makes the schemes C04DS1 and C04DS2 outperform C04 and C04DS3. Figures 2-3 show the corresponding computations using the sixth-order and eight-order accurate schemes. These results are qualitatively similar to the fourth order accurate results in Fig. 1. As in the fourth-order case, the unsplit centered schemes C06 and C08 and the C06DS3 and C08DS3 split schemes are insensitive to the choice of source term, and show the longest stability for three of the source terms. However, with the source term, e G , the split schemes DS1,DS2, and DS4 all outperform DS3 and the centered schemes, also for the higher-order computations. The blow-up of the solution occurs at a somewhat earlier time with the higherorder accuracy. The conclusions here are only tested on one particular 2D smooth flow test case. It is possible that different performance would be obtained for multiscale smooth flows, flows with discontinuities, shock-free turbulent flows and turbulence with strong shocks.
Furthermore, only four choices of the Ducros et al. like splitting among many choices were chosen. Although it would not be hard to define other different split schemes, but there is no unique guiding principle for how this should be done to maximize the stability.

MHD test case with shocks -2D Orzag-Tang vortex test case
The Orzag-Tang vortex starts from initial data ρ = 25/9 (u, v, w) = (− sin y, sin x, 0) p = 5/3, and is solved in two space dimensions on a domain of size 2π × 2π with periodic boundary conditions. The solution will be computed up to time 3.14, when the smooth initial data will have developed shock waves. Because shock waves appear, this problem requires nonlinear shock capturing filter as post processing to the non-dissipative Ducros et al. split form of the base scheme. Results will be displayed as density contour levels at time 3.14 together with contours of div B at the same time. All contour plots in this section use the same values of the contour levels. Furthermore, a logscale plot of the evolution of the norm of div B will also be given. The computational domain was discretized by a grid with 100 × 100 grid points. Figures 4-6 show results by the standard sixth-order centered base scheme filtered with the dissipation of the fifth-order WENO scheme after each time step (C06+WENO5fi), for the source term choices e = e 0 , e = e G , and e = e J respectively.
The centered base scheme preserves the discretized div B perfectly as indicated on the middle subplots of Figs. 4-6 provided the nonlinear filter was not applied to the three magnetic field equations of the MHD. The empty middle subplots in Figs. 4-6 indicate an error-free discrete div B to the order of the truncation error based on contour levels used to make the plot. It is noted that for uniform grid, when the base scheme solved for the full MHD equations (include all three magnetic field equations) and the nonlinear filter step is not applied on the three magnetic field equations, we do not need the standard DivB cleaning procedure. Due to a page limitation, see on expanded version of this paper for the error contours (not empty middle subplots) when   the filter step were employed on all three magnetic fields [22]. We also would like to mention that the entropy conserving schemes developed recently by Sjögreen & Yee [23,24] under the same framework of the nonlinear filter approach might not have the perfect div B preservation property for the same test case for certain choice of MHD governing equation set and entropy numerical fluxes. Warning: For magnetic field dominating flows, switching off the nonlinear filter on the magnetic field components might not be stable for problems where the shock waves are very strong. In this case, good standard div(B) cleaning method is needed.    In these computations, the nonlinear filter was not applied to the magnetic field components, giving perfect div B preservation.
The Orzag-Tang problem does not involve long time integration, and therefore the long time stability studied in the Alfvén wave problem, is not an issue here. Any advantage of the split schemes over unsplit schemes deriving from the long time integration properties, would not be visible in this problem. The Orzag-Tang problem is solved here to show that the split base schemes are performing well also on problems where discontinuities are present.
DNS computations of 2D and 3D MHD turbulence to illustrate the stability and accuracy of our proposed schemes are a subject of our current work in progress.

Conclusions
The Yee & Sjögreen and Kotov et al. [33,15,14] high order numerical method with the Ducros et al. skew-symmetric type of splitting for compressible gas dynamics is extended to ideal MHD. Four Ducros et al. type skew-symmetric splittings together with four formulations of the MHD are proposed and investigated. The proposed methods are able to stabilize non-dissipative high order central finite-difference methods for hyperbolic conservation laws. These methods are also able to maintain their high accuracy and stability for rapidly developing flows and also can improve the stability and accuracy of long time wave propagations and long time integration of turbulent flows. It was shown that the choice of the div B source term in the MHD equations has a large influence on the nonlinear stability. However, for problems where shock waves are present, without added dissipation the computed solutions exhibit oscillations. The nonlinear filter method by Yee & Sjögreen was demonstrated to maintain stability and give highly accurate computed solutions when using skew-symmetric splitting as a pre-processing step for the central base scheme. The support of the DOE/SciDAC SAP grant DE-AI02-06ER25796 is acknowledged. The work was performed with the third author as a research scientist at the Bay Area Envirornmental Research Institute. Financial support from the NASA TTT/RCA program for the second author is gratefully acknowledged. The authors are grateful to Dr. Alan Wray of NASA Ames Research Center for the numerous valuable discussion throughout the course of this work.