- About this Journal ·
- Abstracting and Indexing ·
- Aims and Scope ·
- Annual Issues ·
- Article Processing Charges ·
- Articles in Press ·
- Author Guidelines ·
- Bibliographic Information ·
- Citations to this Journal ·
- Contact Information ·
- Editorial Board ·
- Editorial Workflow ·
- Free eTOC Alerts ·
- Publication Ethics ·
- Reviewers Acknowledgment ·
- Submit a Manuscript ·
- Subscription Information ·
- Table of Contents
The Scientific World Journal
Volume 2013 (2013), Article ID 672187, 14 pages
Nonoscillatory Central Schemes for Hyperbolic Systems of Conservation Laws in Three-Space Dimensions
Department of Mechanical Engineering, The University of Akron, Akron, OH 44325-3903, USA
Received 22 April 2013; Accepted 18 June 2013
Academic Editors: Z. Guo and S. A. Mohiuddine
Copyright © 2013 Andrew N. Guarendi and Abhilash J. Chandy. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
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.
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
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–3]. The first-order Lax-Friedrichs scheme  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 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), , which was further extended to higher orders of accuracy [5–7], as well as to multidimensional systems [8–10]. The main ingredient in the construction of the NT method was a second-order nonoscillatory, monotonic upstream scheme for conservation laws-(MUSCL-) type , 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 . However, applying the fully discrete NT scheme to the second-order convection-diffusion equations did not provide the desired resolution of discontinuities (see [12–14]) due to the accumulation of excessive numerical dissipation [12, 15]. This led Kurganov and Tadmor  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  and central weighted essentially nonoscillatory (CWENO) reconstruction  (also see [16, 17] for essentially nonoscillatory (ENO) implementations). Such weighted essentially nonoscillatory (WENO) reconstructions were introduced first in an upwind framework  after which they were extended to a central framework [5, 8, 9, 19]. More recently Balbás and Tadmor [20, 21] presented extensions of the semidiscrete central schemes of  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  for systems of equations.
Starting with a general hyperbolic conservation law in three-space dimensions, (1), let the sliding averages of over the cells (see Figure 1) at time level, , be where , and are the cell widths in the -, -, and -directions, respectively.
Following , 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 , , , , and 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 .
The cell interface values in the -, -, and -directions, shown in the previous equations, are defined as
They are calculated via a nonoscillatory piecewise polynomial reconstruction given by the polynomials are determined so that satisfies properties such as conservation of cell averages, accuracy, and nonoscillatory behavior (see ). For more details about the derivation of the polynomials and their properties, see . 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 of the solution and it is reprojected back onto the original grid cells.
The resulting semi-discrete scheme in the limit as is as follows:
with numerical fluxes
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.
3. 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 . 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 , , and . In the -direction they are as follows: where the linear polynomials
conserve the pair of cell averages , and , , respectively, and the parabola centered around ,
is determined so as to satisfy where
Note that for the central parabola equation (10) additional corrections in the and directions guarantee third-order accuracy when the point value is recovered from the neighboring cell averages.
The conservation of the cell averages and the accuracy property (property in [8, 13]) are guaranteed [8, 21] by any symmetric choice of weights (e.g., , ). The nonoscillatory behavior (property in [8, 21]) is attained with nonlinear weights prevents the denominator from being zero (), and the smoothness indicators provide a local measure of the derivatives , switching automatically to the second-order reconstructions and 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 component of and
represents its norm over the discretized solution domain.
Similar reconstructions to (8) can be carried out in the - and -directions and the interface values , and 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.
4. 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 in all the cases presented here.
4.1. 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 equation is solved in a computational domain of size with a total of up to points, and the initial conditions are such that the variable is given by
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 , , , , and . The solution with resolution is highly inaccurate due to excessive dissipation and can be disregarded in this case. As we go from to , the solution starts to improve, and eventually, it does not change; that is, the difference between the and 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 , 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 solution of (19) using the polynomial reconstructions without diagonal smoothing, respectively, for 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 and diffuses as well, as time progresses.
4.2. 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 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, , and the parallel efficiency, , where , , 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 , while the other is with 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 simulation can be run is , while the corresponding number for the simulation is . 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 and processors for the and simulations, respectively.
Figure 7 shows the simulation time for time steps on a log scale, where the point corresponding to a single processor was in fact extrapolated from the nearest point assuming ideal efficiency (). The CPU time decreases linearly with the number of processors which is encouraging. On the same figure speedup and efficiency are close to ideal (red dashed line), with efficiency values ranging between and .
4.3. Euler System of Equations: Richtmyer-Meshkov Instability
Euler equations of gas dynamics are given by
Here and are scalar quantities representing the mass density and total internal energy, respectively. is the velocity field with Euclidean norm . The pressure, , is coupled to the total internal energy, .
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 experiment of Collins and Jacobs  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 grid points on a domain of . For test purposes and in order to have a higher resolution, the domain size here in these simulations is more than smaller in the -direction as compared to experiments.
The initial conditions were adapted from the Mach air/ experimental shock tube configuration of Collins and Jacobs . The adiabatic exponent corresponding to an air mixture was used. The ratio of densities is given by . The initial sinusoidal interface had preshock amplitude cm and wavelength cm. An initial diffusion layer thickness of was used , where the thickness function is if , if and otherwise. , and ( is machine zero). Figure 8 shows the initial condition in terms of density for this case. See  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 cm (maximum in the streamwise direction) and at the ghost points, which is exact and does not generate spurious noise .
Figures 9 and 10 show the instantaneous contour slices of density and isosurfaces of density, respectively, at times given by , , , and . As the RMI instability develops, spikes of 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, .
4.4. 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 , and and represent the magnetic field and its norm, respectively. The pressure, , is coupled to the total internal energy, . Furthermore, the system of MHD equations is augmented by the solenoidal constraint; that is, if the condition is satisfied initially at , then by (24) it remains invariant in time.
The problem, that is investigated here, is a MHD problem, the Orszag-Tang-type problem . 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 , where . Again, grid independence is demonstrated in Figure 11 through density contours on slices across the centerlines planes on coarse () and fine grids (). The results presented here in this section are those computed on the grid using the CWENO reconstruction without diagonal smoothing.
A way of demonstrating the accuracy of a numerical method is to determine whether the solenoidal constraint is maintained throughout the simulation. Since 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 is not equal to ). The schemes presented here when first introduced in  in and 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 on all the -surfaces at a certain instant in time. The maximum value of the magnitude of at this instant is , 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 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 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.
Extensions of the semi-discrete schemes of Balbás and Tadmor  to 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 3rd-order 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 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.
- E. Tadmor, “Non-oscillatory central differencing for hyperbolic conservation laws,” Journal of Computational Physics, vol. 87, no. 2, pp. 408–463, 1990.
- R. J. LeVeque, Numerical Methods for Conservation Laws, Lectures in Mathematics, Birkhäuser, Basel, Switzerland, 1992.
- M. E. Hussaini, B. E. VanLeer, and J. E. VanRosendale, Upwind and High-Resolution Schemes, Springer, New York, NY, USA, 1997.
- K. O. Friedrichs and P. D. Lax, “Systems of conservation equations with a convex extension,” Proceedings of the National Academy of Sciences of the United States of America, vol. 68, pp. 1686–1688, 1971.
- D. Levy, G. Puppo, and G. Russo, “Central WENO schemes for hyperbolic systems of conservation laws,” Mathematical Modelling and Numerical Analysis, vol. 33, no. 3, pp. 547–571, 1999.
- X.-D. Liu and E. Tadmor, “Third order nonoscillatory central scheme for hyperbolic conservation laws,” Numerische Mathematik, vol. 79, no. 3, pp. 397–425, 1998.
- D. Levy, G. Puppo, and G. Russo, “A fourth-order central weno scheme for multidimensional hyperbolic systems of conservation laws,” SIAM Journal on Scientific Computing, vol. 24, no. 2, pp. 480–506, 2003.
- D. Levy, G. Puppo, and G. Russo, “Compact central WENO schemes for multidimensional conservation laws,” SIAM Journal on Scientific Computing, vol. 22, no. 2, pp. 656–672, 2001.
- D. Levy, G. Puppo, and G. Russo, “Third order central WENO scheme for 2D conservation laws,” Applied Numerical Mathematics, vol. 33, no. 1, pp. 415–421, 2000.
- G.-S. Jiang and E. Tadmor, “Nonoscillatory central schemes for multidimensional hyperbolic conservation laws,” SIAM Journal on Scientific Computing, vol. 19, no. 6, pp. 1892–1917, 1998.
- B. VanLeer, “Towards the ultimate conservative difference scheme. V. A second-order sequel to Godunov's method,” Journal of Computational Physics, vol. 32, no. 1, pp. 101–136, 1979.
- A. Kurganov and E. Tadmor, “New high-resolution central schemes for nonlinear conservation laws and convection-diffusion equations,” Communications on Pure and Applied Mathematics, vol. 50, pp. 753–771, 1997.
- A. Kurganov and P. Rosenau, “Effects of a saturating dissipation in Burgers-type equations,” Communications on Pure and Applied Mathematics, vol. 50, no. 8, pp. 753–771, 1997.
- A. Kurganov, D. Levy, and P. Rosenau, “On Burgers-type equations with nonmonotonic dissipative fluxes,” Communications on Pure and Applied Mathematics, vol. 51, no. 5, pp. 443–473, 1998.
- A. Kurganov and D. Levy, “A third-order semidiscrete central scheme for conservation laws and convection-diffusion equations,” SIAM Journal on Scientific Computing, vol. 22, no. 4, pp. 1461–1488, 2001.
- A. Harten, B. Engquist, S. Osher, and S. R. Chakravarthy, “Uniformly high order accurate essentially non-oscillatory schemes, III,” Journal of Computational Physics, vol. 71, no. 2, pp. 231–303, 1987.
- W. Shu, “Numerical experiments on the accuracy of ENO and modified ENO schemes,” Journal of Scientific Computing, vol. 5, no. 2, pp. 127–149, 1990.
- G.-S. Jiang and C.-W. Shu, “Efficient implementation of weighted ENO schemes,” Journal of Computational Physics, vol. 126, no. 1, pp. 202–228, 1996.
- L. Cai, J.-H. Feng, and W.-X. Xie, “A CWENO-type central-upwind scheme for ideal MHD equations,” Applied Mathematics and Computation, vol. 168, no. 1, pp. 600–612, 2005.
- J. Balbás, E. Tadmor, and C.-C. Wu, “Non-oscillatory central schemes for one- and two-dimensional MHD equations: I,” Journal of Computational Physics, vol. 201, no. 1, pp. 261–285, 2004.
- J. Balbás and E. Tadmor, “Nonoscillatory central schemes for one- and two-dimensional magnetohydrodynamics equations. II; High-order semidiscrete schemes,” SIAM Journal on Scientific Computing, vol. 28, no. 2, pp. 533–560, 2006.
- A. Kurganov and G. Petrova, “A third-order semi-discrete genuinely multidimensional central scheme for hyperbolic conservation laws and related problems,” Numerische Mathematik, vol. 88, no. 4, pp. 683–729, 2001.
- C. W. Shu, “Total-variation-diminishing time discretizations,” SIAM Journal on Scientific Computing, vol. 9, pp. 1073–1084, 1988.
- S. Gottlieb, C.-W. Shu, and E. Tadmor, “Strong stability-preserving high-order time discretization methods,” SIAM Review, vol. 43, no. 1, pp. 89–112, 2001.
- R. Richtmeyer, “Taylor instability in shock acceleration of compressible fluids, Comm,” Pure and Applied Mathematics, vol. 13, pp. 297–319, 1960.
- E. E. Meshkov, “Instability of the interface of two gases accelerated by a shock wave,” Fluid Dynamics, vol. 4, no. 5, pp. 101–104, 1972.
- B. D. Collins and J. W. Jacobs, “PLIF flow visualization and measurements of the Richtmyer-Meshkov instability of an air/SF6 interface,” Journal of Fluid Mechanics, vol. 464, pp. 113–136, 2002.
- M. Latini, O. Schilling, and W. S. Don, “Effects of WENO flux reconstruction order and spatial resolution on reshocked two-dimensional Richtmyer-Meshkov instability,” Journal of Computational Physics, vol. 221, no. 2, pp. 805–836, 2007.
- S. A. Orszag and C.-M. Tang, “Small-scale structure of two-dimensional magnetohydrodynamic turbulence,” Journal of Fluid Mechanics, vol. 90, no. 1, pp. 129–143, 1979.
- R. B. Dahlburg and J. M. Picone, “Evolution of the Orszag-Tang vortex system in a compressible medium. I. Initial average subsonic flow,” Physics of Fluids B, vol. 1, no. 11, pp. 2153–2165, 1989.
- G.-S. Jiang and C.-C. Wu, “A high-order WENO finite difference scheme for the equations of ideal magnetohydrodynamics,” Journal of Computational Physics, vol. 150, no. 2, pp. 561–594, 1999.
- P. Arminjon and R. Touma, “Central finite volume methods with constrained transport divergence treatment for ideal MHD,” Journal of Computational Physics, vol. 204, no. 2, pp. 737–759, 2005.
- R. Touma and P. Arminjon, “Central finite volume schemes with constrained transport divergence treatment for three-dimensional ideal MHD,” Journal of Computational Physics, vol. 212, no. 2, pp. 617–636, 2006.