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

Abstract

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.

1. 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–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–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 [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–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 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.

2. Formulation

Starting with a general hyperbolic conservation law in three-space dimensions, (1), let the sliding averages of 𝑢 over the cells 𝐶_{𝑖,𝑗,𝑘}∶=[𝑥_𝑖−1/2,𝑥_𝑖+1/2]×[𝑦_𝑗−1/2,𝑦_𝑗+1/2]×[𝑧_𝑘−1/2,𝑧_𝑘+1/2] (see Figure 1) at time level, 𝑛, be 𝑢^𝑛_{𝑖,𝑗,𝑘}1∶=Δ𝑥Δ𝑦Δ𝑧_{𝐶𝑖,𝑗,𝑘}𝑢𝑥,𝑦,𝑧,𝑡^𝑛𝑑𝑥𝑑𝑦𝑑𝑧, (2) where Δ𝑥,Δ𝑦, and Δ𝑧 are the cell widths in the 𝑥-, 𝑦-, and 𝑧-directions, respectively.

Figure 1 

Modified central differencing in three-space dimensions.

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 𝑎^𝑛_{𝑖+1/2,𝑗,𝑘}∶=max𝜌𝜕𝑓𝑢𝜕𝑢^{LC𝑖+1,𝑗,𝑘},𝜌𝜕𝑓𝑢𝜕𝑢^{RC𝑖,𝑗,𝑘},𝑏^𝑛_{𝑖,𝑗+1/2,𝑘}∶=max𝜌𝜕𝑔𝑢𝜕𝑢^{BC𝑖,𝑗+1,𝑘},𝜌𝜕𝑔𝑢𝜕𝑢^{TC𝑖,𝑗,𝑘},𝑐^𝑛_{𝑖,𝑗,𝑘+1/2}∶=max𝜌𝜕ℎ𝑢𝜕𝑢^{BaC𝑖,𝑗,𝑘+1},𝜌𝜕ℎ𝑢𝜕𝑢^{FC𝑖,𝑗,𝑘}, (3) 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].

Figure 2 

Reconstructions in the 𝑥-, 𝑦-, and 𝑧-directions.

The cell interface values in the 𝑥-, 𝑦-, and 𝑧-directions, shown in the previous equations, are defined as 𝑢^RC_{𝑖,𝑗,𝑘}∶=𝑝^𝑛_{𝑖,𝑗,𝑘}𝑥_𝑖+1/2,𝑦_𝑗,𝑧_𝑘,𝑢^LC_{𝑖,𝑗,𝑘}∶=𝑝^𝑛_{𝑖,𝑗,𝑘}𝑥_𝑖−1/2,𝑦_𝑗,𝑧_𝑘,𝑢^TC_{𝑖,𝑗,𝑘}∶=𝑝^𝑛_{𝑖,𝑗,𝑘}𝑥_𝑖,𝑦_𝑗+1/2,𝑧_𝑘,𝑢^BC_{𝑖,𝑗,𝑘}∶=𝑝^𝑛_{𝑖,𝑗,𝑘}𝑥_𝑖,𝑦_𝑗−1/2,𝑧_𝑘,𝑢^FC_{𝑖,𝑗,𝑘}∶=𝑝^𝑛_{𝑖,𝑗,𝑘}𝑥_𝑖,𝑦_𝑗,𝑧_𝑘+1/2,𝑢^BaC_{𝑖,𝑗,𝑘}∶=𝑝^𝑛_{𝑖,𝑗,𝑘}𝑥_𝑖,𝑦_𝑗,𝑧_𝑘−1/2. (4)

They are calculated via a nonoscillatory piecewise polynomial reconstruction given by 𝑅𝑥,𝑦,𝑧;𝑢^𝑛=_{𝑖,𝑗,𝑘}𝑝^𝑛_{𝑖,𝑗,𝑘}𝑥,𝑦,𝑧𝟏_{𝐶𝑖,𝑗,𝑘}; (5) 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 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: 𝑑𝑑𝑡𝑢_{𝑖,𝑗,𝑘}(𝑡)𝐻=−^{𝑥𝑖+1/2,𝑗,𝑘}(𝑡)−𝐻^𝑥_{𝑖−1/2,𝑗,𝑘}(𝑡)−𝐻Δ𝑥_{𝑦𝑖,𝑗+1/2,𝑘}(𝑡)−𝐻^𝑦_{𝑖,𝑗−1/2,𝑘}(𝑡)−𝐻Δ𝑦_{𝑧𝑖,𝑗,𝑘+1/2}(𝑡)−𝐻^𝑧_{𝑖,𝑗,𝑘−1/2}(𝑡)Δ𝑧 (6)

with numerical fluxes 𝐻^𝑥_{𝑖+1/2,𝑗,𝑘}(𝑡)=12𝑓𝑢^LC_{𝑖+1,𝑗,𝑘}+𝑓𝑢^RC_{𝑖,𝑗,𝑘}−𝑎_{𝑖+1/2,𝑗,𝑘}(𝑡)2𝑢^LC_{𝑖+1,𝑗,𝑘}−𝑢^RC_{𝑖,𝑗,𝑘},𝐻^𝑦_{𝑖,𝑗+1/2,𝑘}(𝑡)=12𝑔𝑢^TC_{𝑖,𝑗+1,𝑘}+𝑔𝑢^BC_{𝑖,𝑗,𝑘}−𝑏_{𝑖,𝑗+1/2,𝑘}(𝑡)2𝑢^TC_{𝑖,𝑗+1,𝑘}−𝑢^BC_{𝑖,𝑗,𝑘},𝐻^𝑧_{𝑖,𝑗,𝑘+1/2}(𝑡)=12ℎ𝑢^FC_{𝑖,𝑗,𝑘+1}+ℎ𝑢^BaC_{𝑖,𝑗,𝑘}−𝑐_{𝑖,𝑗,𝑘+1/2}(𝑡)2𝑢^FC_{𝑖,𝑗,𝑘+1}−𝑢^BaC_{𝑖,𝑗,𝑘} (7)

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 [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 𝑃₋₁(𝑥,𝑦,𝑧), 𝑃₀(𝑥,𝑦,𝑧), and 𝑃₁(𝑥,𝑦,𝑧). In the 𝑥-direction they are as follows: 𝑃^𝑛_{𝑖,𝑗,𝑘}𝑥,𝑦,𝑧=𝑤₋₁𝑃₋₁𝑥,𝑦,𝑧+𝑤₀𝑃₀𝑥,𝑦,𝑧+𝑤₁𝑃₁𝑥,𝑦,𝑧,_{𝑚∈−1,0,1}𝑤_𝑚=1, (8) where the linear polynomials 𝑃₋₁𝑥,𝑦,𝑧=𝑢^𝑛_{𝑖,𝑗,𝑘}+𝑢^𝑛_{𝑖,𝑗,𝑘}−𝑢^𝑛_{𝑖−1,𝑗,𝑘}Δ𝑥𝑥−𝑥_𝑖,𝑃₁𝑥,𝑦,𝑧=𝑢^𝑛_{𝑖,𝑗,𝑘}+𝑢^𝑛_{𝑖+1,𝑗,𝑘}−𝑢^𝑛_{𝑖,𝑗,𝑘}Δ𝑥𝑥−𝑥_𝑖 (9)

conserve the pair of cell averages 𝑢^𝑛_{𝑖−1,𝑗,𝑘}, 𝑢^𝑛_{𝑖,𝑗,𝑘} and 𝑢^𝑛_{𝑖,𝑗,𝑘}, 𝑢^𝑛_{𝑖+1,𝑗,𝑘}, respectively, and the parabola centered around 𝑥_𝑖, 𝑃₀𝑥,𝑦,𝑧=𝑢^𝑛_{𝑖,𝑗,𝑘}−112𝑢^𝑛_{𝑖−1,𝑗,𝑘}−2𝑢^𝑛_{𝑖,𝑗,𝑘}+𝑢^𝑛_{𝑖+1,𝑗,𝑘}+𝑢^𝑛_{𝑖,𝑗,𝑘−1}−2𝑢^𝑛_{𝑖,𝑗,𝑘}+𝑢^𝑛_{𝑖,𝑗,𝑘+1}+𝑢^𝑛_{𝑖+1,𝑗,𝑘}−𝑢^𝑛_{𝑖−1,𝑗,𝑘}2Δ𝑥𝑥−𝑥_𝑖+𝑢^𝑛_{𝑖+1,𝑗,𝑘}−𝑢^𝑛_{𝑖,𝑗,𝑘}+𝑢^𝑛_{𝑖−1,𝑗,𝑘}Δ𝑥²𝑥−𝑥_𝑖², (10)

is determined so as to satisfy 𝑐₋₁𝑃₋₁𝑥,𝑦,𝑧+𝑐₁𝑃₁+𝑥,𝑦,𝑧1−𝑐₋₁−𝑐₁𝑃₀𝑥,𝑦,𝑧=𝑢^𝑛_{𝑖,𝑗,𝑘}+𝑢_{𝑥∣𝑖,𝑗,𝑘}𝑥−𝑥_𝑖+12𝑢_{𝑥𝑥∣𝑖,𝑗,𝑘}𝑥−𝑥_𝑖², (11) where 𝑢^𝑛_{𝑖,𝑗,𝑘}∶=𝑢^𝑛_{𝑖,𝑗,𝑘}−124𝑢^𝑛_{𝑖−1,𝑗,𝑘}−2𝑢^𝑛_{𝑖,𝑗,𝑘}+𝑢^𝑛_{𝑖+1,𝑗,𝑘}+𝑢^𝑛_{𝑖,𝑗−1,𝑘}−2𝑢^𝑛_{𝑖,𝑗,𝑘}+𝑢^𝑛_{𝑖,𝑗+1,𝑘}+𝑢^𝑛_{𝑖,𝑗,𝑘−1}−2𝑢^𝑛_{𝑖,𝑗,𝑘}+𝑢^𝑛_{𝑖,𝑗,𝑘+1},𝑢_{𝑥∣𝑖,𝑗,𝑘}∶=𝑢^𝑛_{𝑖+1,𝑗,𝑘}−𝑢^𝑛_{𝑖−1,𝑗,𝑘},𝑢2Δ𝑥_{𝑥𝑥∣𝑖,𝑗,𝑘}∶=𝑢^𝑛_{𝑖+1,𝑗,𝑘}−𝑢^𝑛_{𝑖,𝑗,𝑘}+𝑢^𝑛_{𝑖−1,𝑗,𝑘}Δ𝑥². (12)

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 𝑢^𝑛_{𝑖−1,𝑗,𝑘},𝑢^𝑛_{𝑖,𝑗,𝑘},𝑢^𝑛_{𝑖+1,𝑗,𝑘} and the accuracy property (property 𝐏_𝟐 in [8, 13]) are guaranteed [8, 21] by any symmetric choice of weights 𝑐_𝑚 (e.g., 𝑐₋₁=𝑐₁=1/4, 𝑐₀=1/2). The nonoscillatory behavior (property 𝐏_𝟑 in [8, 21]) is attained with nonlinear weights 𝑤_𝑚=𝛼_𝑚∑_𝑙𝛼_𝑙,with𝛼_𝑚=𝑐_𝑚𝜖+𝐼𝑆_𝑚²,𝑚,𝑙∈{−1,0,1}. (13)𝜖≪1 prevents the denominator from being zero (𝜖=10⁻⁶), 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 𝐼𝑆₋₁=𝑢^𝑛_{𝑖,𝑗,𝑘}−𝑢^𝑛_{𝑖−1,𝑗,𝑘}²,𝐼𝑆₁=𝑢^𝑛_{𝑖+1,𝑗,𝑘}−𝑢^𝑛_{𝑖,𝑗,𝑘}²,𝐼𝑆₀=133𝑢^𝑛_{𝑖+1,𝑗,𝑘}−2𝑢^𝑛_{𝑖,𝑗,𝑘}+𝑢^𝑛_{𝑖−1,𝑗,𝑘}²,+14𝑢^𝑛_{𝑖+1,𝑗,𝑘}−𝑢^𝑛_{𝑖−1,𝑗,𝑘}². (14)

In the case of systems of equations, the smoothness indicators are given by the norm-scaled average of the componentwise indicators, 𝐼𝑆_𝑚^(𝑛𝑞), given by 𝐼𝑆_𝑚=1𝐸𝑄^𝐸𝑄_𝑛𝑞=11𝑢^(𝑛𝑞)₂+𝜖𝐼𝑆_𝑚^(𝑛𝑞),𝑚∈{−1,0,1}, (15) where 𝑢^(𝑛𝑞) stands for the (𝑛𝑞)th component of 𝑢 and 𝑢^(𝑛𝑞)₂=_{𝑖,𝑗,𝑘}𝑢^(𝑛𝑞)_{𝑖,𝑗,𝑘}²Δ𝑥Δ𝑦Δ𝑧 (16)

represents its 𝑙₂ norm over the discretized solution domain.

The interface values can now be calculated from (8)–(10) as follows: 𝑢^RC_{𝑖,𝑗,𝑘}∶=𝑝^𝑛_{𝑖+1/2,𝑗,𝑘}𝑥,𝑦,𝑧=𝑤₋₁𝑢^𝑛_{𝑖,𝑗,𝑘}+12𝑢^𝑛_{𝑖,𝑗,𝑘}−𝑢^𝑛_{𝑖−1,𝑗,𝑘}+𝑤₀𝑢^𝑛_{𝑖,𝑗,𝑘}−112𝑢^𝑛_{𝑖−1,𝑗,𝑘}−2𝑢^𝑛_{𝑖,𝑗,𝑘}+𝑢^𝑛_{𝑖+1,𝑗,𝑘}+𝑢^𝑛_{𝑖,𝑗−1,𝑘}−2𝑢^𝑛_{𝑖,𝑗,𝑘}+𝑢^𝑛_{𝑖,𝑗+1,𝑘}+𝑢^𝑛_{𝑖,𝑗,𝑘−1}−2𝑢^𝑛_{𝑖,𝑗,𝑘}+𝑢^𝑛_{𝑖,𝑗,𝑘+1}+14𝑢^𝑛_{𝑖+1,𝑗,𝑘}−𝑢^𝑛_{𝑖−1,𝑗,𝑘}+14𝑢^𝑛_{𝑖+1,𝑗,𝑘}−2𝑢^𝑛_{𝑖,𝑗,𝑘}+𝑢^𝑛_{𝑖−1,𝑗,𝑘}+𝑤₁𝑢^𝑛_{𝑖,𝑗,𝑘}+12𝑢^𝑛_{𝑖+1,𝑗,𝑘}−𝑢^𝑛_{𝑖,𝑗,𝑘},𝑢^LC_{𝑖,𝑗,𝑘}∶=𝑝^𝑛_{𝑖−1/2,𝑗,𝑘}𝑥,𝑦,𝑧=𝑤₋₁𝑢^𝑛_{𝑖,𝑗,𝑘}−12𝑢^𝑛_{𝑖,𝑗,𝑘}−𝑢^𝑛_{𝑖−1,𝑗,𝑘}+𝑤₀𝑢^𝑛_{𝑖,𝑗,𝑘}−112𝑢^𝑛_{𝑖−1,𝑗,𝑘}−2𝑢^𝑛_{𝑖,𝑗,𝑘}+𝑢^𝑛_{𝑖+1,𝑗,𝑘}+𝑢^𝑛_{𝑖,𝑗−1,𝑘}−2𝑢^𝑛_{𝑖,𝑗,𝑘}+𝑢^𝑛_{𝑖,𝑗+1,𝑘}+𝑢^𝑛_{𝑖,𝑗,𝑘−1}−2𝑢^𝑛_{𝑖,𝑗,𝑘}+𝑢^𝑛_{𝑖,𝑗,𝑘+1}−14𝑢^𝑛_{𝑖+1,𝑗,𝑘}−𝑢^𝑛_{𝑖−1,𝑗,𝑘}+14𝑢^𝑛_{𝑖+1,𝑗,𝑘}−2𝑢^𝑛_{𝑖,𝑗,𝑘}+𝑢^𝑛_{𝑖−1,𝑗,𝑘}+𝑤₁𝑢^𝑛_{𝑖,𝑗,𝑘}−12𝑢^𝑛_{𝑖+1,𝑗,𝑘}−𝑢^𝑛_{𝑖,𝑗,𝑘}. (17)

Similar reconstructions to (8) can be carried out in the 𝑦- and 𝑧-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.

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 ×Δ𝑡=CFLmax𝑎_{𝑖,𝑗,𝑘}²(Δ𝑥)²+max𝑏_{𝑖,𝑗,𝑘}²Δ𝑦²+max𝑐_{𝑖,𝑗,𝑘}²(Δ𝑧)²⁻¹, (18) where 𝑎_{𝑖,𝑗,𝑘},𝑏_{𝑖,𝑗,𝑘}, and 𝑐_{𝑖,𝑗,𝑘} are the speeds of propagation (3). The CFL number was chosen to be 0.5 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 (𝑢)_𝑡+∇⋅𝑢²=0. (19)

The equation is solved in a 3D computational domain of size 1×1×1 with a total of up to 100³ points, and the initial conditions are such that the variable 𝑢 is given by 𝑢=0.5;0.1≤𝑥≤0.5,0.1≤𝑦≤0.5,0.1≤𝑧≤0.5, (20)

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³, 60³, 80³, 100³, and 120³. The solution with 40³ resolution is highly inaccurate due to excessive dissipation and can be disregarded in this case. As we go from 60³ to 120³, the solution starts to improve, and eventually, it does not change; that is, the difference between the 100³and 120³ solutions is not significant. This also indicates that the solution becomes grid independent.

Figure 3 

Instantaneous solutions, at 𝑋=0.5 and 𝑌=0.5, of the nonlinear scalar equation (19) using the CWENO polynomial reconstruction at time 𝑡=0.2 for various grid resolutions. Lines: dashed black ∶403, solid black ∶603, solid red ∶803, solid blue ∶1003, and solid green ∶1203.

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 4 

Instantaneous solutions, at 𝑋=0.5 and 𝑌=0.5, of the nonlinear scalar equation (19) at time 𝑡=0.2 for 3rd- and 2nd-order CWENO reconstructions and at various grid resolutions. Solid lines: 3rd order, dashed lines: 2nd order, red: 603, blue: 803, and green: 1003.

Figure 5 presents the entire 3D solution of (19) using the polynomial reconstructions without diagonal smoothing, respectively, for 100³ 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.

(a)

(b)

(c)

Figure 5 

Instantaneous solutions of the nonlinear scalar equation (19) using the CWENO polynomial reconstruction diagonal smoothing at times (a) 𝑡=0, (b) 𝑡=0.4, and (c) 𝑡=0.8. Shown are isosurfaces at two values 0.5 and 0.05, along with slices along 𝑋=0.3 and 𝑍=0.3, which are projected to the end of the domain. Initial condition—𝑢=0.5, for 0.1≤𝑥≤0.5, 0.1≤𝑦≤0.5, and 0.1≤𝑧≤0.5.

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 16 processors.

Figure 6 

Mesh and domain decomposition for 64×64×64 grid with a 4×2×2 processor configuration.

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³, while the other is with 256³ 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³ simulation can be run is 16, while the corresponding number for the 256³ 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³ simulations, respectively.

(a)

(b)

(c)

Figure 7 

Parallel scaling analysis for the solution of the Euler hydrodynamics system (21) on two different problem sizes: (a) CPU time for 10 time steps, (b) speedup, and (c) efficiency; solid line with symbols: 1283; dashed line with symbols: 2563. The red dashed line indicates perfect or ideal values.

Figure 7 shows the simulation time for 10 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 (100%). 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 94% and 100%.

4.3. Euler System of Equations: Richtmyer-Meshkov Instability

Euler equations of gas dynamics are given by 𝜌_𝑡+∇⋅𝜌𝝊=0,𝜌𝝊_𝑡+∇⋅𝜌𝝊𝝊^𝑇(𝐸)+𝑝𝐈=0,_𝑡𝛾+∇⋅1𝛾−1𝑝+2𝜌𝑣²𝝊=0. (21)

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, 𝐸=(1/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.9cm³. 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₆ 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 𝜌_SF6/𝜌_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.5cm was used [28], where the thickness function is 𝑆(𝑥,𝑦,𝑧)=1 if 𝑑≤0, =exp(−𝛼|𝑑|⁸) 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.

Figure 8 

Instantaneous contours of initial density across centerline y-direction (𝑌=0) at 𝑡=0ms.

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 𝑡=1ms, 𝑡=2ms, 𝑡=3ms, and 𝑡=4ms. As the RMI instability develops, spikes of SF₆ 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].

(a)

(b)

(c)

(d)

Figure 9 

Instantaneous contours of density across centerline y-direction (𝑌=0) at (a) 𝑡=1ms, (b) 𝑡=2ms, (c) 𝑡=3ms, and (d) 𝑡=4ms.

(a)

(b)

(c)

(d)

Figure 10 

Instantaneous isosurfaces of density (𝐶=7 and 𝐶=9) on a subdomain of size 4.5×4.5×4.5cm3, (a) 𝑡=1ms, (b) 𝑡=2ms, (c) 𝑡=3ms, and (d) 𝑡=4ms.

4.4. Ideal Magnetohydrodynamic (MHD) Equations: Orszag-Tang Vortex System

The system of equations for ideal magnetohydrodynamics (MHD) is given by 𝜌_t+∇⋅𝜌𝝊=0,𝜌𝝊_t+∇⋅𝜌𝝊𝝊^T+1𝑝+2𝐵²𝐈−𝐁𝐁^𝑇(𝐁)=0,_t()(𝐸)−∇×𝝊×𝐁=0,_t𝛾+∇⋅1𝛾−1𝑝+2𝜌𝜐²𝝊−(𝝊×𝐁)×𝐁=0. (22)(23)(24)(25)

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, 𝐸=(1/2)𝜌𝜐²+𝑝/(𝛾−1)+(1/2)𝐵². Furthermore, the system of MHD equations is augmented by the solenoidal constraint; that is, if the condition ∇⋅𝐁=0 is satisfied initially at 𝑡=0, then by (24) it remains invariant in time.

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: 𝜌=𝛾²,,𝑝=𝛾,𝐮=−sin𝑦sin𝑧,sin𝑥sin𝑧,0𝐁=−sin𝑦sin𝑧,sin2𝑥sin𝑧,sin2𝑥sin𝑦 (26)

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³) and fine grids (256³). The results presented here in this section are those computed on the 256³ grid using the CWENO reconstruction without diagonal smoothing.

(a)

(b)

Figure 11 

Instantaneous contours of density across 𝑥-, 𝑦-, and 𝑧-direction centerlines (=𝜋) at 𝑡=0.2 for two different grid sizes 1283 (a) and 2563 (b).

A way of demonstrating the accuracy of a numerical method is to determine whether the solenoidal constraint ∇⋅𝐁=0 is maintained throughout the simulation. Since ∇⋅𝐁=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 ∇⋅𝐁 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 ∇⋅𝐁 on all the 𝑍-surfaces at a certain instant in time. The maximum value of the magnitude of ∇⋅𝐁 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.

Figure 12 

Instantaneous surface plot of the divergence of the magnetic field (∇⋅𝐁) on all the 𝑍-surfaces at 𝑡=0.5 for the solution of the Orszag-Tang system.

Figure 13 

Isosurface of the density at a value 𝜌𝐶=3.0 at 𝑡=0.2.

(a)

(b)

(c)

Figure 14 

Instantaneous contours of density across 𝑥-, 𝑦-, and 𝑧- direction centerlines (=𝜋) at (a) 𝑡=0.2, (b) 𝑡=0.4, and (c) 𝑡=0.8 (grid: 2563).

(a)

(b)

(c)

Figure 15 

Instantaneous slices across the 𝑥-direction centerline (=𝜋) showing vectors of the 2D magnetic field (𝐵𝑦,𝐵𝑧) colored by the magnitude of the magnetic field at (a) 𝑡=0.2, (b) 𝑡=0.4, and (c) 𝑡=0.8 (grid: 2563).

5. 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 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 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.