Nonoscillatory Central Schemes for Hyperbolic Systems of Conservation Laws in Three-Space Dimensions

We extend a family of high-resolution, semidiscrete central schemes for hyperbolic systems of conservation laws to three-space dimensions. Details of the schemes, their implementation, and properties are presented together with results from several prototypical applications of hyperbolic conservation laws including a nonlinear scalar equation, the Euler equations of gas dynamics, and the ideal magnetohydrodynamic equations. Parallel scaling analysis and grid-independent results including contours and isosurfaces of density and velocity and magnetic field vectors are shown in this study, confirming the ability of these types of solvers to approximate the solutions of hyperbolic equations efficiently and accurately.


Introduction
Over the past couple of decades, much work has gone into the construction, analysis, and implementation of modern numerical algorithms for the approximate solution of systems of nonlinear hyperbolic conservation laws of the form ( ) + ( ) + ( ) + ℎ( ) = 0. (1) Numerical solutions of these equations are of tremendous practical importance as they govern a variety of physical phenomena in natural and engineering applications. In particular, a number of high-resolution schemes have been developed and tested for this purpose [1][2][3]. The first-order Lax-Friedrichs scheme [4] is actually the forerunner for a large class of central high-resolution schemes that have seen much development in recent years. This paper presents the formulation and testing of such 3D central high-resolution schemes.
These central schemes enjoy a major advantage of simplicity over methods like upwind schemes, in that no approximate Riemann solvers are involved in their construction. To provide a brief background, in 1990, Nessayahu and Tadmor introduced a fully discrete second-order nonoscillatory central scheme (NT scheme), [1], which was further extended to higher orders of accuracy [5][6][7], as well as to multidimensional systems [8][9][10]. The main ingredient in the construction of the NT method was a second-order nonoscillatory, monotonic upstream scheme for conservation laws-(MUSCL-) type [11], piecewise linear interpolant (instead of the piecewise constant one employed in the Lax-Friedrichs scheme) in combination with a higher-order solver for the time evolution [1]. However, applying the fully discrete NT scheme to the second-order convection-diffusion equations did not provide the desired resolution of discontinuities (see [12][13][14]) due to the accumulation of excessive numerical dissipation [12,15]. This led Kurganov and Tadmor [12] to the development of a set of second-order semi-discrete central schemes, which had smaller dissipation than the original NT scheme, and, unlike the fully discrete central schemes, they could be efficiently used with time steps as small as required by the CFL stability restriction due to the diffusion term.
The schemes presented in Kurganov and Tadmor's study saw further developments in the form of third-order extensions [15] and central weighted essentially nonoscillatory (CWENO) reconstruction [8] (also see [16,17] for essentially nonoscillatory (ENO) implementations). Such weighted essentially nonoscillatory (WENO) reconstructions were introduced first in an upwind framework [18] after which they were extended to a central framework [5,8,9,19]. More 2 The Scientific World Journal recently Balbás and Tadmor [20,21] presented extensions of the semidiscrete central schemes of [12] with arbitrary order, , specifically third-and fourth-order reconstructions with the possibility of additional reconstructions in the diagonal directions. In this paper, we present and test the third-order accurate semidiscrete central schemes of Balbás and Tadmor [21] for 3D systems of equations.

Formulation
Starting with a general hyperbolic conservation law in threespace dimensions, (1), let the sliding averages of over the cells , , Figure 1) at time level, , be where Δ , Δ , and Δ are the cell widths in the -, -, and -directions, respectively. Following [21], the first step passing from a fullydiscrete to a semi-discrete formulation consists of reducing the size of the staggered cells covering the Riemann fans. So cells with variable width of order (Δ ) are used, by incorporating the maximum local speeds of propagation across the cell interfaces [12,22], which are given by where the superscripts LC, RC, BC, TC, BaC and FC stand for the left, right, bottom, top, back, and front centers, respectively ( Figure 2). The previous equations are exact for genuinely nonlinear and linearly degenerate fields [21].
The cell interface values in the -, -, and -directions, shown in the previous equations, are defined as the polynomials , , ( , , ) are determined so that ( , , ; ) satisfies properties such as conservation of cell averages, accuracy, and nonoscillatory behavior (see [21]). For more details about the derivation of the polynomials and their properties, see [21]. This information allows the separation of regions of smoothness and regions of nonsmoothness, where the discontinuities propagate in one or more directions. In the nonsmooth regions, staggered evolution over the respective control volumes is used to obtain cell averages, while, for the smooth regions, the polynomial , , ( , , ) and the corresponding fluxes are integrated over the respective control volume to obtain cell averages in these regions. This way, a new polynomial, say ( , , ; ), is formed, which is reconstructed from these smooth and non-smooth portions The Scientific World Journal   of the solution and it is reprojected back onto the original grid cells.
The resulting semi-discrete scheme in the limit as Δ → 0 is as follows: for the third-order CWENO reconstruction without any diagonal smoothing (diagonal smoothing described in Section 3). There is also an option of carrying out diagonal smoothing, the equations of which are not provided here.

Implementation of Semidiscrete Central Schemes: Cell Interface Reconstructions
This section provides a third-order nonoscillatory reconstruction in three-space dimensions, that was implemented for computing the solutions of hyperbolic conservation laws (see (19), (21), and (22)- (25)). The reconstruction of the point values of presented in this section is the third-order CWENO polynomial reconstruction of Kurganov and Levy [15]. The properties of the piecewise quadratic polynomial, that is used here, were presented in [15,21]. In each cell , , (defined before), the polynomials [ , , ( , , )] in (4) are written as a convex combination of three polynomials −1 ( , , ), 0 ( , , ), and 1 ( , , ). In the -direction they are as follows: , , ( , , ) = −1 −1 ( , , ) + 0 0 ( , , ) where the linear polynomials conserve the pair of cell averages −1, , , , , and , , , +1, , , respectively, and the parabola centered around , is determined so as to satisfy Note that for the central parabola equation (10)  guaranteed [8,21] by any symmetric choice of weights (e.g., −1 = 1 = 1/4, 0 = 1/2). The nonoscillatory behavior (property P 3 in [8,21]) is attained with nonlinear weights ≪ 1 prevents the denominator from being zero ( = 10 −6 ), and the smoothness indicators provide a local measure of the derivatives ( , , ), switching automatically to the second-order reconstructions −1 and 1 in the presence of steep gradients and avoiding the onset of spurious oscillations [9,15,18], and in this case they read In the case of systems of equations, the smoothness indicators are given by the norm-scaled average of the componentwise indicators, ( ) , given by where ( ) stands for the ( )th component of and represents its 2 norm over the discretized solution domain.
The interface values can now be calculated from (8) Similar reconstructions to (8) can be carried out in theand -directions and the interface values FC , , , BaC , , , TC , , , and BC , , can be derived in a straightforward manner. Section 4 presents specific test cases related to the governing equations presented before along with a parallel scaling analysis of the implementation of such schemes on multiple platforms.

Numerical Test Cases
This section presents results from the solutions of a single scalar equation (19), the Euler equations of gas dynamics (21), and the ideal equations of magnetohydrodynamics (22) and systems of equations and in applications of a scalar convection problem, the Richtmyer-Meshkov instability, and the Orszag-Tang vortex problem, respectively. For all the calculations presented here, we choose a uniform grid in physical space. For temporal discretization, the third-order strong stability-preserving (SSP) Runge-Kutta [23,24] is used and the time step is dynamically calculated to satisfy the CFL restriction given by where , , , , , , and , , are the speeds of propagation (3). The CFL number was chosen to be 0.5 in all the cases presented here.

Single Scalar Equation: Inviscid Burgers Equation.
The presentation of three-dimensional semi-discrete schemes, formulated in Section 3, begins with the solution of a single scalar equation given by The Scientific World Journal and zero elsewhere. Firstly, the solution of the previous equation is presented at various grids to assess the convergence of such a scheme. Also note that no diagonal smoothing is applied in this case. In Figure 3, solutions are obtained at various grid resolutions such as 40 3 , 60 3 , 80 3 , 100 3 , and 120 3 . The solution with 40 3 resolution is highly inaccurate due to excessive dissipation and can be disregarded in this case. As we go from 60 3 to 120 3 , the solution starts to improve, and eventually, it does not change; that is, the difference between the 100 3 and 120 3 solutions is not significant. This also indicates that the solution becomes grid independent.
Next, the schemes presented in this paper are compared to the second-order scheme of Kurganov and Tadmor [12], where a nonoscillatory reconstruction consisting of a minmod limiter was used. Figure 4 shows the 3rd-and 2nd-order solutions at three different grid resolutions in order to assess if there is in fact an improvement in resolution with the order of the scheme. With regard to resolving the discontinuities, the 3rd-order scheme presented here exhibits slightly more dissipation than the 2nd-order method at all resolutions. The solution is expected to dissipate more with the application of diagonal smoothing. Figure 5 presents the entire 3D solution of (19) using the polynomial reconstructions without diagonal smoothing, respectively, for 100 3 grid resolution. Isosurfaces at two different values and the slices in two different directions are shown at various times. The nonzero region or "cube" moves towards the corner (1, 1, 1) and diffuses as well, as time progresses.

Parallel Performance Analysis for System of Equations.
The computational requirements for the solution of hyperbolic problems could become prohibitive in the case of three-dimensional, geometrically complex enclosures. These requirements increase further when realistic fluid flows like viscous or turbulent flows are considered, thereby requiring larger computational effort and memory. Recent developments in high-performance computing promise a substantial increase in computational speed and offer new possibilities for more accurate simulations. Three-dimensional domain decomposition is used to speed the calculations, where the computational domain is decomposed into a number of rectangular blocks with each processor being responsible for a single block. An example of this decomposition can be seen by the gaps in the grid in Figure 6 for the specific case of 16 processors.
Most of the calculations in the interior of each of the subdomains are independent of the domain decomposition and can continue as if they are performed serially. Problems arise near the subdomain boundaries where, for example, finite differences calculated adjacent to the subdomain boundaries may need several points outside the subdomain. To support these circumstances, two rows of "ghost points" are carried along with the interior solutions that contain copies of the interior solution from the neighboring subdomain. These points are exchanged and updated from neighboring processors as needed to ensure that all near-wall calculations are performed with current variable values.
If a uniform grid is used, then the subdomains in each direction will contain equal number of grid points. However, for a nonuniform grid, the the division locations between the subdomains need to be selected to provide good load balancing or an equivalent amount of work for each processor in each time step. Hence, for the purpose of a scaling analysis, Figure 7 illustrates the CPU time, parallel speedup, = serial / parallel , and the parallel efficiency, = / , where serial , parallel , and are the CPU time for serial and parallel runs and the number of processors, respectively. The scaling analysis is presented for the numerical solutions of the Euler hydrodynamic system presented before; (21) diagnostics for two different problem sizes are presented: one is with 128 3 , while the other is with 256 3 number of points. The simulations were conducted on the Glenn cluster at Ohio Supercomputing. Due to the high memory requirements of the code, the lowest number of processors on which a 128 3 simulation can be run is 16, while the corresponding number for the 256 3 simulation is 32. In order to present a complete scaling analysis, that is, to calculate speedup and efficiency, it is assumed that these quantities are ideal up to 16 and 32 processors for the 128 and 256 3 simulations, respectively.
Here and are scalar quantities representing the mass density and total internal energy, respectively. = ( , V, ) is the velocity field with Euclidean norm 2 := ‖ ‖ 2 . The pressure, , is coupled to the total internal energy, = (1/2) 2 + /( − 1). The evolution of the Richtmyer-Meshkov instability (RMI) [25,26] is considered in this section. RMI arises when a shock passes through an interface between two fluids of widely ranging densities. A generic feature of these systems, as is the case for fluid turbulence in general, is the existence of fluctuations on multiple length scales. Three-dimensional simulations of the reshocked RMI modeled after the Mach 1.21 experiment of Collins and Jacobs [27] are presented in the present work. The simulations use the 3rd-order CWENO reconstruction method without diagonal smoothing (to avoid excess dissipation resulting from it) using 1024× 512 × 512 grid points on a domain of 17.8 × 8.9 × 8.9 cm 3 . For test purposes and in order to have a higher resolution, the domain size here in these simulations is more than 50% smaller in the -direction as compared to experiments. The initial conditions were adapted from the Mach 1.21 air/SF 6 experimental shock tube configuration of Collins and Jacobs [27]. The adiabatic exponent = 1.24815 corresponding to an air mixture was used. The ratio of densities is given by SF 6 / air = 4.063. The initial sinusoidal interface ( , ) = sin(2 / ) sin(2 / ) had preshock amplitude = 0.2 cm and wavelength = 5.933 cm. An initial diffusion layer thickness of = 0.5 cm was used [28], where the thickness function is ( , , ) = 1 if ≤ 0, = exp(− | | 8 ) if 0 < < 1 and 0 otherwise. = ( + ( , ) + − )/(2 ), and = − ln ( is machine zero). Figure 8 shows the initial condition in terms of density for this case. See [28] for a 2D implementation of these initial conditions.
The following boundary conditions were used: (a) inflow at the test section entrance in the streamwise -direction; (b) reflecting at the end wall of the test section in the streamwise direction; and (c) being periodic in the -and -directions corresponding to the cross-section of the test section. The reflecting boundary condition is implemented by reversing the normal component of the velocity vector: ( , ) = − ( , ) at = 17.8 cm (maximum in the streamwise direction) and at the ghost points, which is exact and does not generate spurious noise [28]. Figures 9 and 10 show the instantaneous contour slices of density and isosurfaces of density, respectively, at times given by = 1 ms, = 2 ms, = 3 ms, and = 4 ms. As the RMI instability develops, spikes of SF 6 fall into the air. Following this initial growth, the spikes roll-up and additional complex structures begin to appear. The results presented here are qualitatively similar to those of other studies, for example, [28].

Ideal Magnetohydrodynamic (MHD) Equations: Orszag-
Tang Vortex System. The system of equations for ideal magnetohydrodynamics (MHD) is given by Here and are scalar quantities representing the mass density and total internal energy, respectively, = ( , , ) is the velocity field with Euclidean norm 2 := ‖ ‖ 2 , and B = ( 1 , 2 , 3 ) and 2 := ‖B‖ 2 represent the magnetic field and its norm, respectively. The pressure, , is coupled to the total internal energy, = (1/2) 2 + /( − 1) + (1/2) 2 . Furthermore, the system of MHD equations is augmented by the solenoidal constraint; that is, if the condition ∇ ⋅ B = 0 is satisfied initially at = 0, then by (24)  The 3D problem, that is investigated here, is a 3D MHD problem, the Orszag-Tang-type problem [29]. The evolution of the vortex system involves the interaction between several shock waves traveling at various speed regimes [30,31], which makes the problem especially attractive for numerical experiments. The initial data for this problem are the following: with 0 ≤ , , ≤ 2 , where = 5/3. Again, grid independence is demonstrated in Figure 11 through density contours on slices across the centerlines planes on coarse (128 3 ) and fine grids (256 3 ). The results presented here in this section are those computed on the 256 3 grid using the CWENO reconstruction without diagonal smoothing.
The Scientific World Journal 13 A way of demonstrating the accuracy of a numerical method is to determine whether the solenoidal constraint ∇ ⋅ B = 0 is maintained throughout the simulation. Since ∇ ⋅ B = 0 initially, theoretically it should remain so throughout the simulation. However, the accumulation of numerical errors can usually lead to nonphysical phenomenon known as magnetic monopoles (when ∇ ⋅ B is not equal to 0). The schemes presented here when first introduced in [21] in 1D and 2D frameworks did not require an explicit enforcement of the solenoidal constraint for producing stable and reasonably accurate solutions, and hence no such treatment is used here either. Figure 12 shows the surface plots of ∇ ⋅ B on all thesurfaces at a certain instant in time. The maximum value of the magnitude of ∇ ⋅ B at this instant is 0.41, which is actually representative of the entire simulation. Although these simulations did not show any kind of instabilities, the divergence values are reasonably large for such computations and go to show that in 3D some divergence treatment [32,33] might be necessary. Figure 13 shows a density isosurface, and Figure 14 shows contours of density on three slices across the = = = planes. Figure 15 shows the 2D magnetic field vector colored by magnetic field magnitude. These results are similar to those of previous studies [32,33] and demonstrate the ability of such higher-order central schemes to resolve the shocks that the vortex system develops while maintaining the simplicity and ease of implementation typical of this black-box type of finite difference schemes.

Conclusions
Extensions of the semi-discrete schemes of Balbás and Tadmor [21] to 3D are presented and tested for the first time in this paper. The numerical test cases chosen include evolution of (a) a single scalar equation or an inviscid Burgers equation, (b) the Richtmeyer-Meshkov instability (RMI) using Euler hydrodynamic equations, and (c) the Orszag-Tang vortex system using ideal magnetohydrodynamic equations. Grid independence was demonstrated for two of the three cases presented here. For the single scalar equation test case, our 3rdorder results exhibited more dissipation than those with corresponding 2nd-order schemes with minmod reconstructions. Parallel scaling analysis showed almost ideal efficiencies and speedups based on the assumption that they hold ideal values up until 16 processors. The results obtained with these schemes for the Richtmeyer-Meshkov instability and the Orszag-Tang vortex system confirm the ability of this type of solver to approximate the discontinuous solutions of Eulerian gas dynamics and ideal magnetohydrodynamics equations.