Numerical Dissipation Control in High-Order Methods for Compressible Turbulence: Recent Development ^†

Abstract

This comprehensive overview presents our continued efforts in high-order finite difference method (FDM) development for adaptive numerical dissipation control in the long-time integration of direct numerical simulation (DNS), large eddy simulation (LES), and implicit LES (ILES) computations of compressible turbulence for gas dynamics and MHD. The focus is on turbulence with shock wave numerical simulations using the adaptive blending of high-order structure-preserving non-dissipative methods (classical central, Padé (compact), and dispersion relation-preserving (DRP)) with high-order shock-capturing methods in such a way that high-order shock-capturing methods are active only in the vicinity of shock/shear waves, and high-gradient and spurious high-frequency oscillation regions guided via flow sensors. Any efficient and high-resolution high-order shock-capturing methods are good candidates for the blending of methods procedure. Typically, the adaptive blending of more than one method falls under two camps: hybrid methods and nonlinear filter methods. They are applicable to unstructured finite volume, finite element, discontinuous Galerkin, and spectral element methods. This work represents the culmination of over 20 years of high-order FDM developments and hands-on experience by the authors and collaborators in adaptive numerical dissipation control using the “high order nonlinear filter approach”. Extensions of these FDM versions to curvilinear nonuniform, freestream-preserving moving grids and time-varying deforming grids were also developed. By examining the construction of these two approaches using the high-order multistage type of temporal discretization, the nonlinear filter approach is made more efficient and less CPU-intensive while obtaining similar accuracy. A representative variety of test cases that compare the various blending of high-order methods with standalone standard methods is illustrated. Due to the fact that our nonlinear filter methods are not well known in compressible turbulence with shock waves, the intent of this comprehensive overview is for general audiences who are not familiar with our nonlinear filter methods. For readers interested in the implementation of our methods into their computer code, it is hoped that the long overview will be helpful.

1. Preliminary, Specifics, Objectives, and Outlines

Compressible turbulent flows consist of a disparity of space and time scales during different stages of turbulence development. These flows are usually characterized by very large Reynolds numbers and even larger magnetic Reynolds numbers for magnetized turbulence in magnetohydrodynamics (MHD). DNS, LES, and ILES computations of compressible turbulence require time accurate nonlinearly stable methods suitable for long-time integration with minimal use of numerical dissipation. Designing methods with a minimal amount of numerical dissipation while maintaining numerical stability and accuracy is an intricate balancing act. For compressible turbulence with shock waves, contrary to rapidly developing unsteady flows, standard high-order shock-capturing methods are either too dissipative or encounter nonlinear instability. In order to minimize nonlinear instability and aliasing errors while maintaining high accuracy for different evolutions of flow characteristics, DNS, LES, and ILES simulations require new paradigms to construct the desired efficient methods.

For combustion and reacting flows containing stiff source terms and discontinuities, it is a well-known phenomenon that a wrong propagation speed for discontinuities is obtained when conservative numerical methods that were developed for flows without source terms are used [1,2,3,4,5,6,7,8,9,10,11]. The construction of numerical methods for a stable and accurate simulation of turbulence with shear and shock (viscous shock) waves, and to obtain a correct propagation speed for discontinuities in the presence of stiff source terms, shares one important requirement—the minimization of numerical dissipation while maintaining numerical accuracy and stability. The dual requirements to achieve both numerical stability and minimal numerical dissipation are often conflicting since the majority of past and present shock-capturing methods were designed mainly to be robust in flows without turbulence and/or stiff source terms.

To improve the numerical stability, accuracy, and efficient simulation of compressible turbulence for gas dynamics and MHD, we need methods that (a) are structure-preserving or physics-preserving (this includes but it is not limited to numerically preserving the positivity of pressure and density, the preservation of the divergence of the magnetic field for MHD simulations, entropy conservation, momentum conservation, and kinetic energy preservation), (b) can numerically handle stiff source terms correctly if present, (c) can maintain/improve numerical stability while minimizing the added numerical dissipation for the long-time integration of flows containing both shock-free turbulence regions and turbulence with shocks/shear regions during different stages of turbulence development, and (d) are efficient and the most suited to high-performance modern supercomputers. This work represents an outgrowth of the authors’ involvement in developing adaptive numerical dissipation control in high-order methods for compressible turbulence for gas dynamics and MHD with hands-on experience since the late 1990s [12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30].

Even by taking advantage of recent high-order numerical method developments, using one single method during the entire flow evolution has usually resulted in shortcomings in complex turbulence flow development. This is due to the fact that standalone recent/past methods do not possess the aforementioned four properties. A blending of more than one method to cater to the particular flow regions or different stages of turbulence development, guided by a smart flow sensor at each time step on the type and amount of numerical dissipation, is needed, and it would help improve the overall stability, accuracy, efficency, and reliability of simulations. Typically, the adaptive blending of more than one method falls under two camps: hybrid methods and nonlinear filter methods.

Remark 1.

When using the same non-dissipative, high-order spatial discretization blending with the same high-order shock-capturing method, and a flow sensor in standard turbulence test cases, it was shown by authors [14,15] that the nonlinear filter approach is more efficient and less CPU-intensive while obtaining similar accuracy compared to the hybrid method used in the Johnsen et al. [31] report. See examples shown in later sections.

A more recent trend in turbulence/shock numerical simulations uses the adaptive blending of a high-order skew-symmetric form of structure-preserving non-dissipative spatial central methods (classical central and Padé (compact) and DRP) with the desired high-order shock-capturing methods guided via a smart flow sensor in a manner to further minimize the aid of added numerical dissipation. See the 2011 work of Pirozzoli [32,33] for a review of the adaptive blending of structure-preserving methods (kinetic-preserving) and shock-capturing methods for high-speed flows. Less-known earlier work of Yee et al., Vinokur & Yee, Sandham et al., Sjögreen & Yee, and Yee & Sjögreen [16,17,24,25,26,28] employs a structure-preserving entropy-splitting method in conjunction with their nonlinear filter approaches. Skew-symmetric forms of arbitrary high-order structure-preserving methods will be discussed in later sections, and details can be found in [18,19,20,21,22,23,26,29,30]. The work focuses on the Yee & Sjögreen, Sjögreen & Yee et al. “ high order nonlinear filter approach” with skew-symmetric forms of structure-preserving methods [14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30]. The formulation of arbitrary high-order structure-preserving non-dissipative central discretizations is included in a curvilinear grid in later sections for readers interested in the implementation of our methods into their computer code. Preliminarily, the specifics to be addressed and relevancy are discussed below.

1.1. Specifics and Relevancy

This work addresses only the method-of-line (MOL) approach—a semi-discrete approach to spatial discretizations. Proper temporal discretizations are important but are beyond the scope of the current overview. It is assumed that proper time discretizations with time step size control to ensure accuracy in time are used. The considered high-order spatial discretizations are applicable, in principle, to mix time and spatial discretizations of the Lax–Wendroff/MacCormack and ADER types of methods [34]. Most of the shock-capturing spatial discretizations were developed for inviscid flows in the sense of a vanishing viscosity approach, as pure inviscid flows are inherently chaotic [35]. The methods discussed are for standard hyperbolic conservation laws of the compressible Euler equations and for the ideal MHD equations. If physical viscosities are present, it is assumed that they are high-Reynolds-number Navier–Stokes equations. Methods for combustion and reacting flows are also briefly reviewed.

Historically, most computational fluid dynamics (CFD) numerical methods’ development starts with finite difference methods (FDMs) and finite volume methods (FVMs). For simplicity and ease of presentation, only FDMs in a curvilinear nonuniform grid are presented [24,25,26,36,37,38]. Extensions to time-varying nonuniform deforming grids have been developed for thermally perfect gas. For high-order freestream-preserving, moving, and deforming grids via classical central, Padé, and DRP methods, see Yee et al. (1999, 2000), Vinokur & Yee (2000), Sjögreen et al. (2014, 2020), and Yee & Sjögreen (2023) [19,23,24,25,26,30,36,39].

Applicability of Current FDM Development to Other Formulations: Some of the hybrid methods discussed here have already been extended to FVM, DG, FEM, and spectral element methods. Aside from FDM, FVM, DG, and spectral element methods, other classes of CFD methods are particle methods, LBM (the Lattice–Boltzmann method), and gas-kinetic BGK methods. Regardless of the types of spatial discretizations, the blending of non-dissipative or low-dissipative discretization with some form of discontinuous capturing methods is needed for more accurate simulations of compressible turbulence, rather than employing a standalone high-order shock-capturing method in the presence of shocks. For unstructured grids, DG has been dominating the CFD method development field. See [40] for an overview of DG development. Although our nonlinear filter FDM approaches are applicable to FVM, DG, and spectral element methods, fewer practitioners use these more efficient approaches that exhibit similar accuracy to the hybrid approaches [14,15,31].

The original “Residual Distribution Methods” in FDM form were originally developed for time-marching to the steady state but were recently extended to time-accurate flows for unstructured grids. See review articles [41,42,43]. The current FDM development employing the blending of more than one method for compressible turbulence is also applicable to unstructured residual distribution methods.

Particle methods in CFD are numerical tools for the solution of equations of fluid dynamics obtained by replacing the fluid continuum with a finite set of particles. One of the key attributes of particle methods is that pure advection is treated exactly. The Boltzmann equation is used to approximate macroscopic fluid phenomena as mesoscopic shock events of fictitious particles. For gas–kinetic BGK methods, see [44] for an overview. The current FDM development for compressible turbulence is also applicable to the gas-kinetic BGK methods.

Spectral Element Methods: For complicated geometries using unstructured grids, aside from unstructured FVM, DG, and FEM, spectral element methods were initially developed for spectral-like accuracy in computing incompressible flows. They are essentially finite-element methods that use high-degree polynomial basis functions with compact support. For examples of such methods used in CFD, see nek5000 [45], which solves incompressible Navier–Stokes equations, and nektar++, which has a lot more models than nek5000. For turbulence with shocks, hybridized methods between the spectral element method and DG have been gaining popularity for the last decade.

“Spectral difference/spectral volume” methods are the FDM/FVM version of the spectral element method. They are in the same spirit as “spectral element” methods by using local basis functions, but one might achieve “spectral-like” accuracy using high-degree polynomials. See [46] for a short overview.

1.2. Recent Development in High-Order Shock-Capturing and Shock-Fitting Methods

More than a hundred variants of essentially non-oscillatory (ENO) nonlinear shock-capturing methods have appeared in the literature over the last three decades. Moreover, at least triple the number of weighted ENO (WENO) improvements have been proposed, including very CPU-intensive weighted compact nonlinear schemes (WCNS) with shortcomings in that they can only be parallelized plane by plane. A small percentage of them offer the correction of the shortcoming of the original ENO and WENO methods, and some of them focus on solving other type of fluid flow-governing equations. However, a large percentage of them are minor improvements over the original methods for the compressible Euler equations for perfect gas, yet with complicated logical operations and an increase in the CPU time. Furthermore, these variants can sometimes be less stable or are still too dissipative for long-time integrations. Regardless of their performance, they can be part of the recent trend in the blending of non-dissipative discretization with these more CPU-intensive variants of shock-capturing methods. It is the users’ choice to incorporate their own favorite methods.

1.3. Alternative Methods for Shock Capturing—Front-Tracking Methods and Moretti’s Shock-Fitting Methods

The blending of shock-fitting methods with adaptive dissipation control high-order methods presented here should be applicable to compressible turbulence simulations. For an overview of front-tracking methods, see [47]. For an overview of Moretti’s shock-fitting methods, see [48].

1.4. Artificial Viscosity Methods

The artificial viscosity methods alter the governing equations by adding artificial physical dissipation terms into the original PDEs with tuning parameters. They are the oldest methods for discontinuity capturing. The recent work of Margolin & Lloyd-Ronning [49] presents a comprehensive overview of and a new perspective on artificial viscosity methods. The method of Cook & Cabot (2005) [50] also belongs to such an approach. The guiding philosophy behind the hyperviscosity method is to first regularize the equations through the addition of diffusive terms based on artificial viscosity properties and then to solve the regularized equations using a high-order accurate scheme. A key property is the discrimination between vortical and dilatational structures through the shear and bulk viscosities, respectively. The amount of artificial shear and bulk viscosity, thermal conductivity, and mass diffusivity is controlled via tunable coefficients. A fixed set of values has to be tuned according to the type of flow physics. Studies in the literature indicate that the shock-capturing capability of the Cook & Cabot method is not very consistent in always capturing shock turbulence interaction accurately when compared with the blending of more than one method. See Johnsen et al. [31] for some illustrations.

1.5. Riemann Solvers of Euler Equations

For systems of hyperbolic conservation laws, in order to account for all characteristic waves, shock-capturing methods are usually applied to each of the local characteristic fields. This approach is denoted in the literature as the Riemann solver approach. In the 1970s and 1980s, the Godunov–Riemann solver was used for shock-capturing methods. Since the Roe approximate Riemann solver was first published in 1981 [51], it has been heavily used due to immediate advocacy by Ami Harten and others. The Roe solver is linearized, for which it requires an entropy fix, but it is complete, which means that it represents all characteristic waves. See Yee et al. [38,52] for the earliest multi-dimensional entropy fix to avoid carbuncle phenomena, especially in grid-aligned instability for blunt body flows. Many simplified versions that do not account for all characteristic fields were proposed. The local Lax–Friedrich flux (LLF), taking into account the maximum of the characteristic wave only, and the Harten, Lax, and van Leer (HLL) approximate Riemann solvers were first proposed by Harten et al. [53] in 1983. They do not take into account all characteristic waves, and they are less CPU-intensive but more diffusive. The HLLC (Harten, Lax, van Leer, and Contact) Riemann solver is a modified version of the HLL solver that accounts for the presence of intermediate waves, such as contact discontinuities and shear waves. It is more accurate than HLL, as it can solve exactly stationary contact discontinuities. See [54,55] for an overview.

Remark 2.

In 2D numerical flux formulations for the compressible Euler equation, certain FDM shock-capturing methods using some form of Riemann solver formulations can be written in an equivalent FVM forms. See, e.g., the early work of Yee [38]. For the implementation in curvilinear grids of these FDM shock-capturing methods using a Roe–Riemann solver, see Yee et al. (1985) and Yee (1989) [37,38].

Other Approaches to Constructing Riemann Solvers and Alternatives: The 1982 Osher–Solomon solver has some limitations. Its complexity, CPU intensity, and applicability are the three major drawbacks. See [54,55] for details. Recently, a numerical version has been proposed by Dumbser & Toro, 2011a [56], and Dumbser & Toro, 2011b [57]. Alternative approaches are the flux vector splitting (FVS) method for the compressible Euler flux derivative using upwind discretizations initiated by Steger & Warming [58]. The FVS upwind method is an alternative to good Riemann solvers, but it is more diffusive. Some later FVS developments were proposed that require less computational effort. See [55] for an overview.

1.6. Objectives and Outline

The objective of this work focuses only on the specifics and applicability of other methods of construction under the umbrella indicated previously. A summary of these method-blending approaches is presented in the next two sections for FDM. Some of the early developments can be found in [24,25,26,38,59,60]. Section 3 describes the rationale for the need for efficient, non-dissipative high-order structure-preserving methods blending with high-order shock-capturing methods for compressible turbulence. Section 4 gives a short overview of numerical methods for problems with source terms and problems containing stiff source terms and discontinuities to minimize incorrect numerical solutions. Section 4 also illustrates some numerical examples for test cases with or without source terms. Section 5 and Section 6 describe the entropy-conserving split method, together with its recent extensions. Section 7 shows numerical results, comparing entropy-split high-order methods with seven other high-order methods. Section 8 describes our recent generalization to a wider class of entropy-split methods for compressible ideal MHD.

Comparisons between selected methods and the selected variety of test cases are included throughout. The test cases were chosen to illustrate the effectiveness of the blending of selected numerical method approaches. The selected high-order non-dissipative spatial methods are classical central, Padé (compact), and dispersion relation property (DRP) approximations. For illustration purposes, we only considered total variation diminishing (TVD), high-order WENO variants of Jiang & Shu and Balsara & Shu [61,62], as shock-capturing methods used for the blending of methods when comparing with the standard, standalone high-order method approaches. Unless indicated, Roe’s approximate Riemann solver was used for all test cases. From the wide variety of studied test cases, we found that the performance of the blending of more than one method is highly dependent on the flow physics. Regardless of the choice of blending of high-order methods, these test cases illustrate the improved stability and accuracy over standard standalone high-order method approaches.

It is emphasized that these current approaches can also improve stability and accuracy in rapidly developed unsteady flows containing steep gradients and shear and shock waves with the benefit in minimizing the use of added numerical dissipations as well.

2. Overview of Recent Trends in High-Order Methods’ Development for Compressible Turbulence

For the last two decades, various new numerical methods have been developed that cater to simulating shock-free turbulence and turbulence with shock and shear waves. The major current trends in high-order method developments for the further improved stability, accuracy, and predictability of compressible turbulent developments are summarized in the next subsections, based on the authors’ opinion and decades of observation at specialized conferences and workshops and in the literature.

2.1. Rationale for the Need of Efficient High-Order Methods with Smart Flow Sensors for Adaptive Numerical Dissipation Control

In the compressible turbulence modeling and simulation community, it is known that even employing high-order shock-capturing methods that were designed for sharper discontinuity capturing with rapidly developing flows either eventually encounters instability or is too dissipative for DNS and LES/ILES compressible turbulent flow computations, especially for turbulence with discontinuities. See earlier work, e.g., [25,26,63,64]. The rationale is briefly discussed below.

The original ENO and WENO methods utilizing the shock-resolving capability of TVD schemes together with higher-order accuracy were developed in the late 1980s and early 1990s. High-order multiresolution wavelet methods are part of high-order accuracy development. In the past two decades, many variants of the ENO and WENO methods were derived, yet they were still too dissipative away from shock waves for the long-time integration of turbulent flows. For this reason, the straightforward application of these original and recently developed high-order variants of the ENO and WENO methods alone appears not to be a preferred choice for highly resolved turbulence simulations that push the resolution to the limit without resorting to using an extremely fine grid. From now on, the terms TVD, ENO, and WENO refer to their original formulations and their variants.

2.1.1. Recent Trends

To further the effort in the minimization of numerical dissipation away from shock and shear waves, newer turbulence/shock simulations use adaptive blending of high-order ENO or WENO methods with high-order non-dissipative methods in such a way that high-order ENO or WENO schemes are active only near the shock wave and high gradient regions guided via flow sensors. The blending of more than one method dates back to around five decades ago using the second-order method blending with first-order methods. The best-known ones are the flux-corrected transport (FCT) method of Boris and Book [59] in the 1970s and the artificial compression methods of Harten [60].

High-resolution shock-capturing methods, including higher than the third order of accuracy, actually in a loose sense use the blending of more than one method in their construction. This is due to the fact that the majority of high-resolution TVD, ENO, and WENO shock-capturing methods can be written as “high-order central discretization and nonlinear shock-capturing terms using one or more flux limiters or flow variable limiters” (limiters applied to a set of independent characteristics, primitive or conservative variables). See Harten and Yee et al. for such an early procedure of construction [37,38,65] in detail. As a matter of fact, the original TVD method of Harten, via construction at the outset, is written in a second-order central discretization and a controlled second-order dissipation term consisting of limiters away from discontinuities and high-gradient regions. Yet this original procedure for the construction of the well-known TVD types method is not really known in the computational fluid dynamics (CFD) community. There were isolated reports in this century that considered this approach as a new idea for developing new shock-capturing methods.

The following discussion of the blending of more than one method precludes the method of constructing high-order shock-capturing methods in the sense of the aforementioned FCT and TVD scheme type of construction, blending central and nonlinear dissipation terms using limiters at the outset. Instead, we took a high-order shock-capturing method as a nonlinear scheme that was then blended with a non-dissipative high-order linear scheme (e.g., central, Padé, or DRP) with the help of flow sensors to guide the blending of methods.

For the last two decades, the blending of more than one method has basically concentrated on two camps of such development: (a) the hybrid method that switches between high-order non-dissipative methods and high-order shock-capturing methods (e.g., high-order WENO or ENO) guided via a flow sensor for the computed data [32] and (b) the high-order nonlinear filter method of Yee et al., Sjögreen & Yee, and Yee & Sjögreen [17,25,27,28,66,67]. For nonlinear filter methods, after the completion of a full time step of the primary high-order non-dissipative linear scheme, a smart local flow sensor would examine the computed data by flagging the locations and the amount of numerical dissipations that are needed. Then, the computed data would be filtered via a dissipative portion of high-order shock-capturing methods, according to the flagged locations and the percentage amount of shock-capturing dissipation. More sophisticated forms with more than one type of numerical dissipation and local sensors are reported in [14,15,27,28]. For the performance of the nonlinear filter methods using high-order non-dissipative Padé spatial methods as the primary (base) methods, see [21,68]. Many of the recent developers of more complex variants of ENO and WENO consider the aforementioned hybridized approach with their own flow sensors.

These two general approaches, if with similar flow sensors, usually provide a similar accuracy. However, the nonlinear filter approach is more efficient. For example, consider the dimension-by-dimension Riemann solver approach. Only ONE evaluation of a dissipative portion of a chosen high-order shock-capturing method is needed per dimension, per time step, and per grid spacing for the nonlinear filter approach, regardless of the time integrator. Contrary to the nonlinear filter method, for the hybrid method, depending on the time integrator, the number of high-order shock-capturing method evaluations per dimension, per time step, and per grid spacing would be four by a fourth-stage Runge–Kutta time integrator. If physical viscosities are present, the amount of numerical dissipation needed would be automatically adjusted accordingly, using the flow sensor via the nonlinear filter at the current time step instead of at the previous time step using the hybrid method approach. In addition, the numerical solution might experience a non-smooth transition at the switch between different types of schemes using the hybrid method. For complex shock wave and shear surface interactions, the switch mechanism can become less trivial, and frequent switching between two types of schemes can further promote numerical instability beyond the induced instability from the inherent strong nonlinearity and the presence of multiscale physical processes that are dominant features of the subject flow in question.

2.1.2. Nomenclature

Before further discussion, some nomenclature must be defined. Here a “nonlinear filter” is different from the standard spectral or high-order linear filter used to filter the computed data after each time step (or after more than several time steps). Nonlinear filters make use of the dissipative portion of a nonlinear spatial discretization. Here, we use the term nonlinear methods to mean spatial discretizations using some sort of limiters or adaptive choices for discretized polynomials to suppress the spurious oscillation of computed data. Higher than first-order shock-capturing methods are nonlinear methods that are at least higher than the first order away from discontinuities. They are nonlinear methods, as the resulting difference equations are nonlinear in the flow variables even when applied to a linear governing equation set. This is opposed to high-order non-dissipative linear methods (and high-order linear filters) for which the resulting difference equations are linear in the flow variable. In particular, for improved stability, we advocate for the use of dissipative portions of fourth-order or higher high-resolution shock-capturing methods as the major components of the nonlinear filter term.

2.1.3. High-Order Non-Dissipative Linear Spatial Discretizations

Commonly used high-order non-dissipative methods include the standard central, Padé, and DRP methods. Except for higher than fourth-order Padé spatial discretization, high-order non-dissipative discretizations are highly efficient, with a low CPU operations count. The Padé spatial finite difference method was first introduced by Kreiss and Oliger and implemented by Hirsh [69,70], and it was then popularized by Lele [71]. Since Padé spatial discretizations are global discretizations, the disadvantage of the Padé finite difference scheme is that it is an implicit form of spatial discretization, and it lacks flexibility due to the complex matrix inversion, especially for multi-dimensional cases on a non-uniform grid. Padé methods are not highly parallelizable. See [68,72,73]. For high-order Padé spatial discretizations with non-periodic physical boundary conditions, more careful numerical handling of the boundary conditions is needed. For standard high-order central spatial discretizations, well-known summation-by-part (SBP) stable boundary operators are readily available for non-periodic physical boundary conditions. DRP methods [74,75] were designed for computational acoustics. DRP methods can be provided with SBP stable boundary operators similar to standard-centered discretizations. See Sjogreen & Yee [18] for a computational acoustics study. From here on, the terms “Padé discretization (method)” and “compact discretization (method)” are used interchangeably and assumed to describe central Padé spatial discretization. A large volume in the literature uses the term “Compact Method”.

The following gives an overview of the blending of hybrid methods and their nonlinear filter counterparts.

2.2. Short Overview of Hybrid Methods

To improve the shock-capturing capability while maintaining high accuracy away from discontinuities and high gradient regions in the flow, various hybrid methods have been developed since the mid-1990s. There are many variants of hybrid methods. Typical hybrid methods switch between a high-order non-dissipative spatial (central/Padé/DRP) method and a chosen shock-capturing method. For the Padé method, this is most often with the help of added linear high-order compact filters for stabilization.

Consider a 1D hyperbolic conservation law where $\mathbf{U}$ is the vector of flow variables, ${\mathbf{F}}_x$ is the inviscid flux derivative, and t and x represent time and the computational space in x in a computational domain.

$${\mathbf{U}}_t+{\mathbf{F}}_x=0$$

Assume the semi-discrete form uses a uniform FDM grid, $x_j=(j-1)\Delta x$, $j=1,\dots ,N$, where $\Delta x$ is the grid spacing. Denote ${\mathbf{U}}_j$ and ${\mathbf{F}}_j$ as the approximation of $\mathbf{U}$ and $\mathbf{F}$ at grid spacing j. A semi-discrete form of a chosen method is

$$\frac{d{\mathbf{U}}_j}{dt}+\frac{{\mathbf{H}}_{j+1/2}-{\mathbf{H}}_{j-1/2}}{\Delta x}=\mathbf{0}$$

(1)

where ${\mathbf{H}}_{j+1/2}$ denotes numerical fluxes consistent with $\mathbf{F}$. A simple hybrid method switching between two methods can be written as

$${\mathbf{H}}_{j+1/2}=(1-{\nu }_{j+1/2}){\mathbf{H}}_{j+1/2}^C+{\nu }_{j+1/2}{\mathbf{H}}_{j+1/2}^S,$$

where ${\mathbf{H}}_{j+1/2}^C$ is a high-order non-dissipative numerical flux, and ${\mathbf{H}}_{j+1/2}^S$ is a high-order shock-capturing numerical flux. ${\nu }_{j+1/2}$ is a local flow sensor. For a 1D system of hyperbolic conservation laws using a complete set of characteristic waves, the high-order shock-capturing methods are preferred to locally discretize each of the characteristic waves using a complete set of Riemann solvers [24,25,26,37,38]. If a multidimensional method is not used, it is also very common to apply a 1D Riemann solver method, dimension by dimension, for a multi-dimensional system of hyperbolic conservation laws [14,15,17,24,27,28,37].

It was shown that, for compressible turbulence without source terms, hybrid methods in general offer more advantages over the more complicated improved high-order shock-capturing methods alone [14,15,17,24,25,26,27,28].

2.3. Padé vs. Central Non-Dissipative Spatial Differencing

High-order Padé (compact) central spatial differencing, together with a high-order linear compact filter standalone (without hybridization using a shock-capturing method), is currently the method of choice for incompressible flows. For compressible turbulence flows, the hybridized methods would switch between the Padé high-order central spatial differencing and high-order shock-capturing methods with the aid of a flow sensor. For turbulence with shock waves, linear high-order compact filters are also needed, as well as blending with a shock-capturing method in order to stabilize the simulations due to Gibbs’s phenomenon [21,68]. Padé spatial discretizations are global methods in the sense that their influence is felt throughout the computational domain, whereas all standard central spatial discretizations are local methods with a local influence. The use of the upwind biased Padé method blended with shock-capturing methods was proposed in the literature. See, for example, the work of [69]. In general, the use of an upwind-biased Padé method is more dissipative and, thus, less accurate but more stable than the central Padé method. Furthermore, other investigators used the hybridization between high-order standard upwind difference methods and shock-capturing methods. While improving the overall numerical stability, the accuracy away from discontinuities and high gradient flow regions is compromised due to the inherited dissipative nature of upwind or upwind biased Padé discretizations.

Even though Padé spatial discretizations, popularized by Lele, are discretizations with many disadvantages over the classical central discretizations, they are heavily used by practitioners for both incompressible and compressible turbulence simulations. From our investigations, for compressible turbulence containing shock waves employing the blending of high-order non-dissipative Padé discretization with high-order shock-capturing methods, there is no gain in accuracy when compared with the use of the same order of classical central discretization. See [21,68] for some gas dynamics illustrations containing shocks.

To show the performance between the eighth-order central Padé vs. eighth-order classical central discretizations, we illustrate the simulation of a well-known Alf$\acute{\mathrm{v}}$en MHD test case. See [23,30] for details of the problem setup and the ideal MHD equations notation. Figure 1 shows the problem setup. Figure 2 shows the maximum error time evolution comparison. The comparison is between Padé vs. classical central and Padé vs. classical central with the entropy-splitting form of the inviscid flux derivative consisting of two different entropy split parameters, $\beta$. “Ce” denotes the eighth-order central, and “Co” denotes the eighth-order Padé. “CeES” denotes the eighth-order central, and “CoES” denotes the eighth-order Padé on the entropy split form of the inviscid flux derivative. The four eighth-order classical central methods can prolong the accuracy much longer than the Padé method counterparts with or without the entropy-splitting form of the discretization [26]. The error stays small for the longest time for CeES at $\beta =2.5$ and for CoES at $\beta =1$. The time discretization is the fourth-order classical Runge–Kutta method. Here, the “entropy splitting form of the inviscid flux derivative” denotes entropy-conserving discretization by employing non-dissipative central discretization in the entropy-splitting form of the inviscid flux derivative. See later sections for a discussion or our original work [19,20,22,23,26,76] for details. This method is denoted as the entropy split method in the later discussion and also in our recent publications.

Remark 3.

High-order central Padé spatial discretizations in conjunction with linear high-order Padé (compact) filters are methods of choice for many incompressible and low-speed turbulent/acoustic flows due to their advantage of requiring a very low number of grid points per wavelength, based on linear model equation analysis. This desired property of high-order central Padé differencing seems to have diminished in both the multiscale systems of compressible gas dynamic and MHD test cases containing shock waves. Without added numerical dissipation and/or an added high-order linear Padé filter, the standalone Padé discretization hybridized with a shock-capturing method would suffer from Gibbs’s phenomenon [21,68]. Moreover, the Padé spatial differencing requires more floating-point operations per time step per grid point, and it is less beneficial for fast parallel computations than the central spatial finite differences. Consequently, the Padé differencing requires added CPU time in a parallel computer framework, especially for multi-dimensional systems [72,73].

2.4. Overview of Nonlinear Filter

The nonlinear filter approach is less known in the aerodynamics CFD application community, as it has not often been described in CFD reference books or textbooks. This section gives a longer overview than the hybrid method approach, including test cases from our extensive work [17,19,20,22,23,27,28,66,67,76].

A nonlinear filter approach consists of a base scheme step and a post-processing nonlinear filter step.

2.4.1. Base Scheme Step

A full time step is advanced using a high-order, non-dissipative, spatially central scheme (classical central, Padé or DRP method). A summation-by-parts (SBP) boundary operator [77,78] and a matching-order conservative high-order free stream metric evaluation for curvilinear grids [24] are used. A high-order temporal discretization, such as the third-order TVD Runge–Kutta or fourth-order Runge–Kutta (RK3 or RK4) method, is used.

2.4.2. Post-Processing (Nonlinear Filter) Step

Instead of the hybrid approaches, to further improve nonlinear stability and accuracy from the non-dissipative spatial base scheme, the computed data are nonlinearly filtered via a dissipative portion of a high-order shock-capturing scheme with a local smart flow sensor. At each grid point, the local flow sensor(s) is employed to analyze the regularity of the computed flow field from the base scheme step. Only strong discontinuity locations would receive the full amount of shock-capturing dissipation. In smooth regions, no dissipation is added unless high-frequency oscillations are detected. In regions with strong turbulence, if needed, a small fraction of the shock-capturing dissipation can be added to improve stability. Note that the nonlinear filter numerical fluxes only involve the inviscid flux derivatives, regardless of whether the flow is viscous or inviscid. If viscous terms are present, for the ease of the SBP boundary closure’s implementation for the viscous flux derivatives, the same inviscid central difference operator used for the first derivative is employed twice for the viscous flux derivatives. For a variety of local flow sensors with the automatic selection of different flow types, see [28].

The nonlinear filter idea was first introduced and tested by Yee et al. [25,26] using an artificial compression method (ACM) by Harten [60] as the flow sensor. Later, smart flow sensors were developed by Sjögreen & Yee, Yee & Sjögreen, and Kotov et al. [14,15,17,27,28]. Unlike the hybrid approach, in the presence of physical viscosity, the primary high-order non-dissipative scheme includes the discretizations on both the inviscid and viscous portions of the governing equation before the nonlinear filtering post-processing step. This is opposed to the hybrid method. It examines the computed flow data from the previous time step information. In addition, switching between two methods that can create instability during frequent switching between two methods is opposed to the nonlinear filter methods that do not switch between methods. The high-order base scheme again can be a standard central, Padé, or DRP method, depending on the flow type in question.

Denote the computed solution with the base scheme step using $U_{j,k,l}^{*}$ of 3D Euler equations for gas dynamics. After the completion of a full time step of the spatial base scheme step at time n, the final update of the solution after the filter step at time $n+1$ is (with the numerical fluxes in the y- and z-directions suppressed, as well as their corresponding y- and z-direction indices on the x inviscid flux suppressed)

$$\begin{array}{c}\hfill U_{j,k,l}^{n+1}=U_{j,k,l}^{*}-\frac{\Delta t}{\Delta x}[H_{j+1/2}^{*}-H_{j-1/2}^{*}],\end{array}$$

(2)

$$\begin{array}{c}\hfill H_{j+1/2}^{*}=R_{j+1/2}{\overline{H}}_{j+1/2}.\end{array}$$

(3)

The nonlinear filter numerical fluxes usually involve the use of field-by-field approximate Riemann solvers. If the Roe type of approximate Riemann solver [51] is employed, for example, the x-filter numerical flux vector $H_{j+1/2}^{*}$ is evaluated at the $U_{j,k,l}^{*}$ solution from the base scheme step. $R_{j+1/2}$ is the matrix of right eigenvectors of the Jacobian of the inviscid flux vector in terms of Roe’s average states based on $U^{*}$. $H_{j+1/2}^{*}$, and $H_{j-1/2}^{*}$ are “filter” numerical fluxes in terms of Roe’s average states based on $U_{j,k,l}^{*}$. Denote the elements of the filter numerical flux vector ${\overline{H}}_{j+1/2}$ with ${\overline{h}}_{j+1/2}^m,m=1,2,...,5$. The element of the filter numerical flux ${\overline{h}}_{j+1/2}^m$ has the form

$${\overline{h}}_{j+1/2}^m=\frac{{\kappa }_{j+1/2}^m}{2}{w}_{j+1/2}^m{\varphi }_{j+1/2}^m.$$

(4)

Here, $w_{j+1/2}^m$ is a flow sensor to activate the nonlinear numerical dissipation portion of a high-order shock-capturing scheme, $\frac{1}{2}{\varphi }_{j+1/2}^m$. The term ${\kappa }_{j+1/2}^m$ represents a locally determined positive parameter that is less than or equal to one, based on the regularity of the computed date. The choice of the parameter $\kappa$ can be different for different flow types, and it is automatically chosen using the local ${\kappa }_{j+1/2}^m$ described in [28]. However, if the computation uses the Ducros et al. flow sensor, we set ${\kappa }_{j+1/2}^m=1$ [79].

Figure 3 shows the nonlinear filter procedure for a system of 1D hyperbolic conservation laws. Figure 4 shows the nonlinear filter procedure using a combination of more than one flow sensor. Figure 4 also shows how to obtain, e.g., the dissipative portion of a seventh-order WENO (WENO7) denoted as ${\varphi }_{j+1/2}^m=g_{j+1/2}^m-b_{j+1/2}^m$. The term $g_{j+1/2}^m$ represents the $mth$ characteristic term of WENO7, and $b_{j+1/2}^m$ represents the $mth$ characteristic term of the eighth-order central method. The nonlinear filter approach using WENO7 is denoted with WENO7fi on some of the later numerical test result presentations.

2.4.3. A Historical Note on Nonlinear Filter Approaches

The original idea of the nonlinear filter approach of Yee et al. (2000) [25] was based on shock-capturing methods that are written into a central discretization portion and a nonlinear shock-capturing portion (dissipative portion). All variants of TVD schemes are already written in this form by design. In order to obtain the dissipative portion of other high-order shock-capturing methods, the dissipative portion of the considered shock-capturing ${\varphi }_{j+1/2}^m=g_{j+1/2}^m-b_{j+1/2}^m$ can be obtained in a similar manner as WENO7fi.

It is noted that the nonlinear filter step described above should not be confused with the LES filtering operation. For an extension of the blending of more than one numerical method to the ideal MHD, see the next section for a brief description, and see the references cited therein for details.

As outlined briefly above, the nonlinear filter schemes are efficient in their construction. The total computational cost for a given error tolerance is significantly lower than for standard shock-capturing schemes or their hybrid cousins of the same order. One important reason for their efficiency is that the nonlinear shock-capturing filter dissipation is applied after each full time step of the base scheme, whereas a standard shock-capturing/hybrid method evaluates the shock-capturing method at each stage, e.g., m stage of the Runge–Kutta (R-K) time-stepping scheme. Hence, the nonlinear filter approach requires only one shock-capturing method evaluation per time step per grid point per dimension, independent of the time discretization involved.

More complicated nonlinear filter approaches are obtained by combining more than one type of numerical dissipation and several different types of flow sensors. See, e.g., [14,15,17,18,27,28].

2.4.4. Accuracy Comparison: Hybrid vs. Nonlinear Filter Approaches

Since the publication of the Johnsen et al. [31] article in the Journal of Computational Physics in 2010, it has been heavily cited. Many researchers considered the results of the hybrid method of switching between two methods to be far superior compared to our nonlinear filter counterpart (Yee & Sjogreen and Sjogreen & Yee [14,17,25,26,27,28]). However, what is indicated by Johnsen et al. [31] is not a fair/correct comparison. If these two approaches to the blending of two methods use the same methods, with the same flow sensor parameter, they should produce a similar accuracy. That version of the hybrid method used a built-in tuning parameter hard-wired into the hybrid code to cater to the different test cases while switching between two methods. For our nonlinear filter method, it was mandated that we use the same input flow sensor scheme parameter for all test cases, as the tuning of the flow sensor parameter was not built in to our computer code. The guideline for the comparison for the 2010 Journal of Computational Physics Johnsen et al. result was to use the same input parameter unless the tuning of the scheme parameter was built in. The true performance of these same test cases can be found in our published work [14,15]. In other words, in order to make a fair comparison between the nonlinear filter method and the hybrid method, it is important to consider the aspect of tunable parameters. These are sometimes used in the design of flow sensors. A fair method comparison should either not allow any tuning of parameters over a set of test problems or allow both methods to be optimally tuned for each test problem. This is sometimes not considered enough in the published literature, leading to only one of the methods being tuned.

The efficiency and accuracy of the nonlinear filter approach for a wide variety of flow problems can be found, e.g., in [12,13,14,15,16,18,19,20,22,27,28,30].

Figure 5 shows the accuracy and CPU comparison between a hybrid method reported in [31] and a nonlinear filter method using the same non-dissipative high-order linear spatial discretization blending with the same shock-capturing method for a 2D shock–vorticity interaction. For this test problem illustrated here, we used the Ducros et al. flow sensor [79].

Figure 5 also shows a comparison among a second-order TVD, a seventh-order WENO (WENO7), a hybrid scheme (switch between eighth-order spatial central scheme and WENO7 using wavelet flow sensor as the switch indicator), and the filter scheme WENO7fi (an eighth-order spatial central base scheme, “D08”, and the dissipative portion of WENO7 using the same wavelet flow sensor to guide where the WENO7 dissipation was applied at the post-processing nonlinear filter step). A second-order Runge–Kutta method was used for the TVD scheme, and the classical fourth-order Runge–Kutta method was used for the rest of the spatial schemes. For this particular simple 2D shock–vorticity interaction test case with a simple weak planar shock without a structure, WENO7, hybrid, and WENO7fi have the same accuracy. However, there is a large gain in CPU time using the nonlinear filter approach for this turbulence-free test case. For turbulence with shocks, there is a more beneficial gain both in accuracy and CPU time for our nonlinear filter schemes over the standard WENO counterparts.

Figure 6 shows a 3D isotropic turbulence with shocklets problem setup. Figure 7 shows the accuracy of our nonlinear filter method for the isotropic turbulence simulation with the smart flow sensor built in. Comparing the resolution using the hybrid method reported by Johnson et al. [31], one can see that they produce the same accuracy. For the same 3D isotropic problem setup as in Johnson et al. [31] and the performance of our nonlinear filter method, see [15].

These two camps of adaptive blending of methods can be used as a starting point to build more appropriate schemes for more challenging turbulent problems in computational physics. It is remarked that it is the users’ choice of their favorite high-order shock-capturing methods for the aforementioned two camps of approaches with DNS and LES/ILES computations. Since most of the numerical method developments might not numerically preserve the desired physical properties, Section 3 below discusses additional methods to help minimize the use of numerical dissipation and numerically preserve some of the physical properties.

3. Rationale for the Need for Efficient High-Order Structure-Preserving Methods Blending with High-Order Shock-Capturing Methods

In classical mechanics, temporal, spatial, or both finite difference discretizations that conserve a certain physical property (or combination of physical properties) of the governing equations are most often referred to, in a broad sense, as structure-preserving numerical methods. See, e.g., [80]. Temporal and spatially finite difference discretizations that discretely conserve fundamental properties of the chosen governing equations are, hereafter, denoted as structure-preserving methods (SPMs). SPMs have been gaining more attention in the long-time integration of turbulent flows. A combination of these SPMs can further increase the reliability, stability, and accurate simulation of fluid flows, especially for turbulent flows. The majority of existing high-order methods, including hybrid or high-order nonlinear filter methods might not discretely preserve the desirable physical-preserving properties (NOT SPM).

In the applied mathematics CFD community, entropy-conserving and entropy-stable methods have been flourishing for the last two decades [19,20,22,23,26,40,81,82]. Symplectic time integrators are used to prevent phase distortion, e.g., Hamiltonian systems or the Korteweg–deVries equations [83]. In the applied mechanics and CFD community, the importance of developing SPMs that conserve one or more fundamental properties of the governing equations (mass, minimizing phase distortion, momentum and physical entropy conserving, the positivity of density and pressure, and kinetic energy preserving, etc.) have been ongoing for the last four decades. See, e.g., [19,20,22,23,26,76,81,84,85,86,87,88,89] for some of the developments. A combination of more than one SPM is also possible. See [21,30,90] for formulations or the next section. Formulations using the DRP and Padé high-order SPM non-dissipative linear methods are presented in [18,21].

3.1. SPM Finite Difference Formulation in Split Form: Skew-Symmetric Splittings of the Compressible Euler Flux Derivatives

When the high-order non-dissipative spatial method discussed previously discretized a class of the skew-symmetric form of the inviscid Euler flux derivative (instead of the unsplit flux derivative), it was shown that the appearance of instabilities or aliasing errors in the long-time integration of turbulence flows can be delayed. The splitting can be achieved for certain flow variables or the entire convective flux derivative system. They are written in a special split form that is equivalent to the unsplit version; i.e., the split and unsplit governing equations are equivalent. However, after finite difference discretizations, certain different split forms can maintain a discrete entropy conservation, momentum conservation, or kinetic energy preservation property or provide a stable L2-like energy norm estimate for smooth solutions. See articles [19,20,22,23,26,76,81,84,85,86,87,89,91] for discussions of the performance of various skew-symmetric splitting approaches in DNS and LES applications. See Sjögreen et al. [20] for the extension of these skew-symmetric splittings to any even-order split approximation formulation for the Euler and MHD equations.

For simplicity or presentation, consider the derivative of the product ${\left(ab\right)}_x$, where a and b are functions of x. A split approximation can be written as

$${\left(ab\right)}_x=\alpha {\left(ab\right)}_x+\gamma a{b}_x+\beta a_{x}b.$$

(5)

before discretization. The parameters $\alpha ,\gamma$, and $\beta$ are chosen to be equivalent to the original ${\left(ab\right)}_x$. A simple split derivative is found by setting $\alpha =\gamma =\beta =1/2$

$${\left(ab\right)}_x=\frac{1}{2}{\left(ab\right)}_x+\frac{1}{2}a{b}_x+\frac{1}{2}{a}_{x}b$$

(6)

The Ducros et al. splitting [87] utilized this particular non-conservative split form in conjunction with the classical central discretization operator on the derivative terms. They then applied mathematical manipulations to the difference operators to make the final discretization into a conservative difference method. Their construction will be discussed shortly.

A split approximation for the derivative of the product of three functions, ${\left(abc\right)}_x$, where a, b, and c are functions of x, can be written as

$$\begin{array}{c}{\left(abc\right)}_x=\alpha {\left(abc\right)}_x+\gamma [a{\left(bc\right)}_x+bc{a}_x]+\beta [b{\left(ac\right)}_x+ac{b}_x]\hfill \\ +\kappa [c{\left(ab\right)}_x+ab{c}_x]+\delta [bc{a}_x+ac{b}_x+ab{c}_x]\hfill \end{array}$$

(7)

before discretization. The parameters $\alpha ,\gamma$, $\beta$, $\kappa$, and $\delta$ were chosen to be still equivalent to the original ${\left(abc\right)}_x$. A well-known split derivative that can be written in conservative form and can improve nonlinear stability is Kennedy & Gruber [91] splitting, which is kinetic energy-preserving for the Euler equation by setting the first four parameters to $1/4$ and $\delta =0$ for the triple product derivative. See [86,89,91] for related developments.

A less-known splitting of the inviscid flux derivative for the nonlinear Euler equation of gas dynamics is the entropy-splitting method of [26,88,92]. The idea of this entropy splitting is that their split formulas can be used to estimate the $L^2$ norm-like or energy norm of the computed solution for periodic and non-periodic boundary conditions. In the recent terminology, it is entropy-stable (in the sense of possessing an energy estimate) for the Euler equations using high-order classical central or DRP central schemes in conjunction with summation-by-parts boundary operators [18,26,88,92,93,94,95]. See [94] for a follow-up study of the subject. The next two subsections give a short overview of the Ducros et al. and entropy-splitting methods.

3.2. Conservative Splitting of Ducros et al. Type

The Ducros et al. split approximation [87] for the Euler equations of gas dynamics starts with the terms of the split form approximated using

$$\frac{1}{2}D\left(ab\right)+\frac{1}{2}D\left(a\right)b+\frac{1}{2}aD\left(b\right),$$

(8)

where D is a centered finite difference operator, and a and b are functions of x. For this split approximation via central, Padé, or DRP discretization, after mathematical manipulations, the resulting method conserves the momentum (see below). However, it is not obvious how to obtain a norm estimate for nonlinear systems.

The key step in the Ducros et al. [87] split approximation is to rewrite (8) in conservation form. For example, employing the second-order operator $D{u}_j=(u_{j+1}-u_{j-1})/(2\Delta x)$, we have

$$\frac{1}{2}D\left(ab\right)+\frac{1}{2}D\left(a\right)b+\frac{1}{2}aD\left(b\right)=\frac{1}{4\Delta x}{\Delta }_{+}\left[(a_j+a_{j-1})(b_j+b_{j-1})\right],$$

(9)

where ${\Delta }_{+}{q}_j=(q_{j+1}-q_j)$.

Moreover, Equation (9) can be generalized to standard-centered difference operators of the $2p^{th}$-order of accuracy:

$$D_p{u}_j=\frac{1}{\Delta x}\sum _{k=1}^p{\alpha }_k^{\left(p\right)}(u_{j+k}-u_{j-k}).$$

(10)

This is due to the fact that the coefficients ${\alpha }_k^{\left(p\right)}$ satisfy

$$\sum _{k=1}^{p}k{\alpha }_k^{\left(p\right)}=\frac{1}{2}\sum _{k=1}^p{\alpha }_k^{\left(p\right)}{k}^{2n+1}=0,n=1,\dots ,p-1.$$

(11)

The right-hand side of the algebraic identity

$$\begin{array}{c}{a}_{j+k}{b}_{j+k}-a_{j-k}{b}_{j-k}+(a_{j+k}-a_{j-k})b_j+a_j(b_{j+k}-b_{j-k})=\hfill \\ (a_{j+k}+a_j)(b_{j+k}+b_j)-(a_j+a_{j-k})(b_j+b_{j-k})\end{array}$$

(12)

is written in a conservative form as

$$\begin{array}{c}(a_{j+k}+a_j)(b_{j+k}+b_j)-(a_j+a_{j-k})(b_j+b_{j-k})=\hfill \\ {\displaystyle \sum _{m=0}^{k-1}}(a_{j-m}+a_{j+k-m})(b_{j-m}+b_{j+k-m})-{\displaystyle \sum _{m=0}^{k-1}}(a_{j-1-m}+a_{j-1+k-m})(b_{j-1-m}+b_{j-1+k-m}).\end{array}$$

(13)

Thus, the conservative form of the split approximation becomes

$$\begin{array}{c}\frac{1}{2}{D}_p\left(ab\right)+\frac{1}{2}{D}_p\left(a\right)b+\frac{1}{2}a{D}_p\left(b\right)=\hfill \\ \frac{1}{\Delta x}\sum _{k=1}^p\frac{1}{2}{\alpha }_k^{\left(p\right)}\left((a_{j+k}{b}_{j+k}-a_{j-k}{b}_{j-k})+a_j(b_{j+k}-b_{j-k})+(a_{j+k}-a_{j-k})b_j\right)\\ =\frac{1}{\Delta x}\sum _{k=1}^p\frac{{\alpha }_k^{\left(p\right)}}{2}\left(\sum _{m=0}^{k-1}(a_{j-m}+a_{j+k-m})(b_{j-m}+b_{j+k-m})\right.\\ \hfill -\left.\sum _{m=0}^{k-1}(a_{j-1-m}+a_{j-1+k-m})(b_{j-1-m}+b_{j-1+k-m})\right)=\frac{1}{\Delta x}(h_{j+1/2}-h_{j-1/2}),\end{array}$$

(14)

where the numerical flux is defined as

$$h_{j+1/2}=\sum _{k=1}^p\frac{1}{2}{\alpha }_k^{\left(p\right)}\sum _{m=0}^{k-1}(a_{j-m}+a_{j+k-m})(b_{j-m}+b_{j+k-m}).$$

(15)

To simplify the formulas of the conservative form of split approximations for systems of equations, define

$${\Theta }_{j+1/2}^{\left(p\right)}(a,b)=\sum _{k=1}^p\frac{1}{2}{\alpha }_k^{\left(p\right)}\sum _{m=0}^{k-1}(a_{j-m}+a_{j+k-m})(b_{j-m}+b_{j+k-m}).$$

(16)

For 3D Euler equations of gas dynamics equations, for simplicity of presentation, consider the x-direction inviscid flux

$$\mathbf{F}={\left[\rho u,\rho u^2+p,\rho uv,\rho uw,(e+p)u\right]}^T,$$

where $\mathbf{F}=\mathbf{F}\left(\mathbf{U}\right)$ is the inviscid flux vector, and $\mathbf{u}={\left[u\right(x,y,z),v(x,y,z),w(x,y,z\left)\right]}^T$ is the velocity vector in the x-, y-, and z-directions. $\rho$ denotes the density, p is the pressure, and e is the total energy. Denote ${\mathbf{u}}_j={(u_j,v_j,w_j)}^T$ as the discretization at the j grid location with the y and z discretization indices suppressed for simplicity. It turns out that 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 obtained via

$${\mathbf{F}}_x{|}_{x=x_j}\approx \left(\begin{array}{c}\frac{1}{2}D{\rho }_j{u}_j+\frac{1}{2}{\rho }_{j}D{u}_j+\frac{1}{2}{u}_{j}D{\rho }_j\\ \frac{1}{2}D{\rho }_j{u}_j^2+\frac{1}{2}{\rho }_j{u}_{j}D{u}_j+\frac{1}{2}{u}_{j}D{\rho }_j{u}_j+D{p}_j\\ \frac{1}{2}D{\rho }_j{u}_j{v}_j+\frac{1}{2}{\rho }_j{v}_{j}D{u}_j+\frac{1}{2}{u}_{j}D{\rho }_j{v}_j\\ \frac{1}{2}D{\rho }_j{u}_j{w}_j+\frac{1}{2}{\rho }_j{w}_{j}D{u}_j+\frac{1}{2}{u}_{j}D{\rho }_j{w}_j\\ \frac{1}{2}D{u}_j(e_j+p_j)+\frac{1}{2}{u}_{j}D(e_j+p_j)+\frac{1}{2}(e_j+p_j)D{u}_j\end{array}\right),$$

(17)

which, via (15), can be written in a conservative form with a numerical flux function:

$${\mathbf{H}}_{j+1/2}=\frac{1}{2}\sum _{k=1}^p{\alpha }_k^{\left(p\right)}\sum _{m=1}^{k-1}\left(\begin{array}{c}({\rho }_{j-m}+{\rho }_{j+k-m})(u_{j-m}+u_{j+k-m})\\ ({\rho }_{j-m}{u}_{j-m}+{\rho }_{j+k-m}{u}_{j+k-m})(u_{j-m}+u_{j+k-m})+p_{j-m}+p_{j+k-m}\\ ({\rho }_{j-m}{v}_{j-m}+{\rho }_{j+k-m}{v}_{j+k-m})(u_{j-m}+u_{j+k-m})\\ ({\rho }_{j-m}{w}_{j-m}+{\rho }_{j+k-m}{w}_{j+k-m})(u_{j-m}+u_{j+k-m})\\ (e_{j-m}+p_{j-m}+e_{j+k-m}+p_{j+k-m})(u_{j-m}+u_{j+k-m})\end{array}\right).$$

(18)

The more compact notation introduced in (16) allows (18) to be rewritten as

$${\mathbf{H}}_{j+1/2}=\left(\begin{array}{c}{\Theta }_{j+1/2}^{\left(p\right)}(\rho ,u)\\ {\Theta }_{j+1/2}^{\left(p\right)}(\rho u,u)+{\Theta }_{j+1/2}^{\left(p\right)}(p,1)\\ {\Theta }_{j+1/2}^{\left(p\right)}(\rho v,u)\\ {\Theta }_{j+1/2}^{\left(p\right)}(\rho w,u)\\ {\Theta }_{j+1/2}^{\left(p\right)}(e+p,u)\end{array}\right).$$

(19)

Unlike the linearized and symmetrized system of Euler equations, where the split approximations lead to a stability estimate, there is no such stability estimate for the $2p^{th}$-order accurate conservative Ducros et al. splitting. Since the Ducros et al. splitting results in a conservative scheme, it is applicable to problems containing discontinuities.

3.3. Special Class of Entropy-Conserving SPM—Entropy Split Method (Entropy-Splitting Approach)

One class of SPMs that the authors have been developing is the entropy split method for the compressible Euler flux derivatives. The entropy split method was originally denoted as the “Entropy Spitting Approach” [24,26]. This is a less-known skew-symmetric splitting of the Euler equation for gas dynamics [26,88,92]. This skew-symmetric splitting made use of Harten’s symmetrizable form of the Euler equations in terms of the entropy variables [81] to obtain a semi-discrete splitting of the Euler equations with discrete entropy stability (in space) via the summation-by-parts approach. The entropy splitting is written in terms of the sum of a conservative portion and a nonconservative portion. If central (classical central, Padé, or DRP) discretizations are used for the interior scheme (interior grid points) and a summation-by-parts boundary scheme (boundary points), the resulting split method is entropy-conservative and entropy-stable. See the recent result of Sjögreen & Yee [18,19,21,22,29,76,94] for the proof and additional extensions. It is important to point out that the Harten [81] and Gerritsen & Olsson entropy-splitting form incorrectly selected the non-physical branch of inequality and was later corrected by Yee et al.; this is 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. [26,88,92] is a splitting that takes a form that is more suitable for the discrete stable energy norm estimate technique, including a boundary scheme estimate for arbitrary order of central spatial schemes. See Yee et al. [26] for the formulation.

The entropy splitting discussed in Harten [81] for the inviscid flux derivative of the 1D Euler equations, $\mathbf{F}{\left(\mathbf{U}\right)}_x$, for a perfect gas occurs via the entropy variables W in the following:

$${\mathbf{F}}_x=\frac{\beta }{\beta +1}{\mathbf{F}}_x+\frac{1}{\beta +1}{\mathbf{F}}_W{W}_x,\beta \ne -1$$

(20)

$$W={[w_1,w_2,w_3,w_4,w_5]}^T=\frac{p^{*}}{p}{[e+\frac{\alpha -1}{\gamma -1}p,-\rho u,-\rho v,-\rho w,\rho ]}^T,$$

(21)

where

$$p^{*}=-{\left(p{\rho }^{-\gamma }\right)}^{\frac{1}{\alpha +\gamma }}$$

(22)

and

$$\beta =\frac{\alpha +\gamma }{1-\gamma },\alpha >0\mathrm{or}\alpha <-\gamma .$$

(23)

See Yee et al. [16,25,26] for the formulation, the choice of $\beta$, and early numerical examples.

In the original entropy-splitting method ( Entropy-Splitting Approach), the Euler flux derivatives approximation is not written in the usual numerical flux form. These entropy split methods are entropy-conserving and entropy-stable, but they are usually not conservative numerical methods without the additional reformulation proposed in, e.g., Sjögreen & Yee [19,20,22,23,96].

See Section 8 and [19,20,22,23,30,96] for the additional development of the compressible Euler and ideal MHD equations.

3.3.1. Key Entropy-Conserving Methods for the Compressible Euler and Ideal MHD Equations

Currently, an entropy-conserving method by Tadmor [82,97] is the method of choice for rapidly developing unsteady flows, especially in the development of unstructured grid methods (e.g., the DG method). Our work on typical test cases using the FDM formulation found that the Tadmor entropy-conserving method [82,97] of the same order requires twice the arithmetic operations as the entropy split methods. In addition, Harten’s entropy function can be part of a Tadmor-type entropy-conserving method family [19,76] with similar accuracy and stability. A comparison among methods can be found in [18,19,76].

3.4. Newer Class of Entropy Split Methods in Numerical Flux Form: Not Relying on the Homogeneous Property of the Inviscid Flux and Symmetrizable Inviscid Flux Derivatives

In 2019, Sjögreen & Yee [22] derived a high-order conservative numerical flux for the non-conservative portion of the entropy splitting of the Euler flux derivatives by taking advantage of the homogeneity property of the inviscid flux and symmetrizable inviscid flux derivatives. Due to the construction, this conservative numerical flux requires more arithmetic operations and is less stable than the original entropy split method (not fully conservative). An alternative simple approach for problems containing shock waves is to use a shock detector to switch the entropy split discretization to a standard centered approximation in the neighborhood of shock waves in order to avoid the wrong shock speed.

Instead of taking advantage of the homogeneity property of the inviscid flux and symmetrizable inviscid flux derivatives, a wider class of entropy split methods that do not require the homogeneous property were developed; see our 2022 paper [23]. The formulation is written in standard numerical flux form. This wider class of the entropy split method was extended for the equations of ideal MHD by the authors in [23]. See Section 8 for a discussion and numerical comparison.

3.5. Methods with Multiple Structure-Preserving Properties

It is our assessment that employing high-order SPM non-dissipative methods in conjunction with the nonlinear filter approach should be a method of choice for high-resolution simulations involving compressible turbulence up to mildly supersonic regimes. For turbulence with shocks, there is an even larger gain both in accuracy and CPU time of the nonlinear filter schemes over their standard WENO counterparts.

Forms of the structure-preserving (SPM) numerical fluxes that belong to the class of skew-symmetric splitting of the inviscid flux derivatives in curvilinear grids are presented in Section 5, Section 6, Section 7 and Section 8. Depending on the type of flow physics, high-order central, DRP, or Padé methods are applied to the skew-symmetric split form of the inviscid flux derivative as the baseline method before the post-processing nonlinear filter step, if needed, especially for problems with shocks. Selected methods with multiple structure-preserving properties are given with numerical examples.

Before the general formulation of SPM numerical fluxes, first, a brief background on numerical methods to avoid/minimize the wrong propagation speed of discontinuities for problems containing stiff source terms and discontinuities developed by Wang et al. and Yee et al. [8,11,98,99] is presented. Some test cases for problems without source terms and problems containing stiff source terms and discontinuities are included in the next section.

4. Methods for Problems with Source Terms and Stiff Source Terms with Shocks

Insufficient spatial/temporal resolution may cause an incorrect propagation speed for discontinuities and nonphysical states using standard dissipative numerical methods that were developed for non-reacting flows. This numerical phenomenon was first observed by Colella et al. [100] in 1986, when they considered both the reactive Euler equations and a simplified system obtained by coupling the inviscid Burgers equation with a single convection/reaction equation. LeVeque and Yee [5] showed that a similar spurious propagation phenomenon can be observed even with scalar equations by properly defining a model problem with a stiff source term. They introduced and studied the simple one-dimensional scalar conservation law with an added nonhomogeneous parameter-dependent source term

$$u_t+u_x=S\left(u\right),$$

(24)

$$S\left(u\right)=-\mu u(u-\frac{1}{2})(u-1).$$

(25)

When the parameter $\mu$ is very large, a wrong propagation speed for the discontinuity phenomenon via dissipative numerical methods will be observed in coarse grid computations. In reacting flows, $\frac{1}{\mu }$ can be described as the reaction time. In order to isolate the problem, LeVeque and Yee solved (24) and (25) via the fractional step method using Strang splitting [101]. For this particular source term, the reaction (ODE) step of the fractional step method can be solved exactly. In their study using a pointwise evaluation of the source term ($S\left(u\right)$ is evaluated at the j grid point index, i.e., $S\left(u_j\right)$ for each time evolution), the phenomenon of a wrong propagation speed for discontinuities is connected with the smearing of the discontinuity caused due to the spatial discretization of the advection term. They found that the propagation error is due to the numerical dissipation contained in the scheme, which smears the discontinuity front and activates the source term in a nonphysical manner. The smearing introduces a nonequilibrium state into the calculation. Thus, as soon as a nonequilibrium value is introduced in this manner, the source term turns on and immediately restores the equilibrium while, at the same time, shifting the discontinuity to a cell boundary. By increasing the spatial resolution by an order of magnitude, they were able to improve the correct propagation speed. It is remarked here that, in a general stiff source term problem, a sufficient spatial resolution is as important as the temporal resolution when the reaction step of the fractional step method cannot be solved exactly.

For the last three decades, the wrong speed phenomenon has attracted a large volume of research work (see, e.g., [1,3,4,6,7,9,102,103,104]). Various strategies have been proposed to overcome this wrong speed difficulty for one to two species cases with a single reaction. Since numerical dissipation that spreads the discontinuity front is the cause of the wrong propagation speed of discontinuities, a natural strategy is to avoid any numerical dissipation in the scheme. In combustion, level set and front-tracking methods were used to track the wavefront to minimize this spurious behavior [1,4,6]. See Wang et al. [8] for a comprehensive overview of the last two decades of development. Wang et al. also proposed a new high-order finite difference method with subcell resolution for advection equations with stiff source terms for a single species to overcome the difficulty. For a subcell resolution method for multi-species, see Wang et al. [105].

Our studies in [11] found that using the Strang form of the operator splitting between the homogeneous part of the system of governing equations (without the source term) and the nonlinear stiff source terms is more stable than solving the fully coupled fluid-reaction equations. Here, Strang operator splitting is different from the skew-symmetric splitting of the inviscid flux derivatives. Strang operator splitting is commonly used in the combustion community to solve reacting flow problems. Figure 8 shows the schematic of the Strang operator splitting. For the subcell resolution method, the subsequent two figures, Figure 9 and Figure 10 show the schematic of subcell resolution methods. For problems with source terms employing any higher than first-order methods, we need to use methods that are well balanced. See [98,99] for a well-balanced WENO method for chemical reacting flows.

In the Strang operator splitting approach, all of the aforementioned numerical methods (structure-preserving methods, shock-capturing methods, and the blending of more than one method) can be used to solve the reacting equations portion without the source terms first. Then, the source term is solved using an ODE (ordinary differential equation) solver, as indicated in Figure 9 and Figure 10.

4.1. Test Cases for Blending of Numerical Methods Using Flows WITH or WITHOUT Source Terms

Here, we select a set of test cases that were performed by the authors showing the performance of the blending of numerical method approaches for flows with or without source terms. These test cases show the blending of numerical methods without applying the skew-symmetric splitting of the inviscid flux derivative. Test cases using different structure-preserving non-dissipative spatial discretizations are discussed as well. It is important to point out that well-balanced shock-capturing methods should be used for problems containing source terms in general. Furthermore, well-balanced shock-capturing methods that can capture the correct propagation speed of discontinuities should be used for problems containing stiff source terms and discontinuities. See some of the well-balanced methods with our collaborators [8,11,98,99,105] that can capture the correct propagation speed of discontinuities.

4.1.1. Chapman–Joguet (C-J) Detonation Containing Stiff Source Terms and Discontinuities

The first example illustrates the effect of numerical dissipation on the propagation speed of a discontinuity in the presence of a stiff source term. We computed a 1D C-J detonation for a one-reaction Arrhenius source term. This example, which was considered in [8], has been extensively studied, e.g., in [3,7,8,11,13]. See [3,7,8,11] for the initial flow conditions and additional results.

Figure 11 shows the 1D problem setup. Figure 12 shows a comparison of four high-order accurate shock-capturing schemes, WENO5, WENO5fi, WENO5fi + split, and our high-order subcell resolution version of WENO5 [8,11] (WENO5/SR). The new WENO5/SR method uses the Strang form of the splitting between the homogeneous part of the system of governing equations (without the source term) and the nonlinear stiff source terms. For each time integration, the first step is to solve the homogenous system of governing equations. Then, a subcell resolution detector is used to correct the wrong numerical speed of propagation of discontinuities before solving the next step of the nonlinear stiff source term. It is remarked that the Strang splitting is not to be confused with the skew-symmetric splitting of the inviscid flux derivative as in the Ducros et al. splitting [87]. The reference solution is computed via the WENO5 scheme on a uniform grid with 10,000 points. Here, WENO5fi + split denotes the use of sixth-order classical central spatial discretization, together with the dissipative portion of WENO5; “split” here denotes the use of the Strang form of the splitting between the homogeneous part of the system of governing equations (without the source term) and the nonlinear stiff source terms in solving the governing equations.

Figure 13 shows a 2D C-J detonation problem setup. Figure 14 shows a 1D cross section of density at $t=1.7\times {10}^{-7}$ via the same four methods as the 1D case on a uniform coarse grid of $200\times 40$. The CFL = $0.05$ and $k_0=0.5825\times {10}^{10}$. The right figure is a close-up of the vicinity of the discontinuity. The reference solution is achieved via WENO5 using $4000\times 800$ uniform grid points. WENO5 is more dissipative and also gives the largest error in the location of the discontinuity. WENO5fi performs better and is less dissipative because the WENO5 dissipation is not used everywhere in the computational domain. WENO5/SR gives the most accurate propagation speed, but WENO5fi + split compares very well with WENO5/SR for this particular problem and grid size. The probable reason is that, because of the stabilizing effect of the splitting, less dissipation is needed to keep WENO5fi + split stable. WENO5/SR and WENO5fi + split capture the correct structure using fewer grid points than in Helzel et al. [3] and in Tosatto & Vigevano [7]. Figure 15 shows the CPU comparison among methods. For a detailed grid refinement, error comparison for different CFL, and stiffness parameter comparisons, see our work [11].

4.1.2. 1D and 2D 13-Species EAST Simulations

The second two examples show the computations of a 1D and 2D 13-species chemical reacting flow by Yee et al. and Kotov et al. [11,13,14].

These 1D and 2D simulations consist of a 13-species chemical reacting flow, which is a simplified problem related to the NASA Electric Arc Shock Tube (EAST) experiment. Figure 16 shows the schematic of the EAST experiment. Figure 17 and Figure 18 show the 1D and 2D problem setup. Figure 19 shows the comparison of three grids using the same WENO5fi method (left) and the comparison of five methods with the reference solution (right). The left subfigure of Figure 19 indicates how the shock/shear locations are dependent on the grid size using the same high-order method. The relative width between the shock and shear is the function of the grid spacing compared with the reference solution. The right subfigure indicates how the shock/shear locations are dependent on five different numerical methods (TVD, TVDafi + split, WENO5-llf, WENO5Pafi + split, and TVD/SR). See [13,14] for the description of the five methods.

Figure 20 shows the temperature and pressure evolution of the 2D 13-species test case. Figure 21 shows the comparison of methods, indicating how the shock/shear locations are dependent on the three methods. Different boundary layers are predicted via the three different methods as well. Here, WENO5-LLF denotes WENO5 using the local Lax–Friedrichs flux (LLF), and WENO5P-LLF denotes a positivity-preserving WENO5 of Zhang & Shu [106] using the LLF. The relative width between the shock and shear is a function of the methods.

4.1.3. Three-Dimensional Supersonic Shock–Turbulence Interaction via a Structure-Preserving Method [14]

The third example is a 3D supersonic turbulence flow simulation using our high-order nonlinear filter method. This is a more challenging test case, and it is very CPU-intensive computation with long transient simulation before obtaining meaningful turbulent statistics. Figure 22 shows the problem setup. To illustrate an instantaneous turbulent pattern, Figure 23 shows the instantaneous velocity field, $u_x$ (top) and $u_y$ (bottom), at a slice where $z=\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{s}\mathrm{t}$, obtained with DNS on a grid of $1553\times {256}^2$ points. A turbulent flow enters over the left boundary at Mach $1.5$. The turbulence interacts with a quasi-steady shock wave near the inflow boundary. The figure shows the x- and y-directions velocity field in a $z=const.$ slice of the three-dimensional computational domain. See [14] for details of the problem set-up and method comparisons. The numerical method is an eighth-order, accurate, non-dissipative method, together with an adaptive nonlinear filter using a dissipation portion of the seventh-order positivity-preserving WENO7 of Zhang and Shu [106] (WENO7Pfi). The Ducros et al. skew-symmetric splitting is the structure-preserving eighth-order central spatial discretization. The turbulent Mach number is 0.16. The incoming turbulence causes wrinkles in the shock front. The classical fourth-order Runge–Kutta method (RK4) is used for temporal discretization. Grid refinement and a comparison of standard high-order shock-capturing methods with the high-order structure-preserving nonlinear filter method were performed in [14].

4.2. Problems with Random Forcing Source Term Simulation

Since all of the aforementioned methods and our nonlinear filter methods were originally developed for turbulent flows without source terms, their straightforward application to problems involving forced turbulence or turbulence with volumetric energy sources can result in spurious numerics [8,9,11,98,99,105]. It was shown in Yee et al. [11] and other works that Strang operator splitting is more stable than solving the full Euler equation with source terms included. For problems with time-varying stochastic source terms, there has not been much numerical method development. Since the forcing term in turbulence simulations with rms Mach numbers up to $0.6$ is not very stiff, spurious numerics associated with the wrong shock speed are very unlikely.

With the above in mind, we used methods developed in [8,11,98,99,105] as a guideline for the simulation of stochastically forced homogeneous turbulence in a periodic domain. The spatially seventh-order nonlinear filter methods with adaptive dissipation control developed by the authors and collaborators [8,11,98,99,105] were employed to solve the homogeneous Euler equations, and an ODE solver was used to solve the stochastic forcing terms via the Strang operator splitting approach [101]. For the majority of numerical tests presented below, time steps based on a Courant–Friedrich–Lewy (CFL) number of $0.8$ were used. The forced turbulent simulations without the Strang operator splitting had to use a much smaller CFL constraint [8,11]. See [107] for some early studies. The results of higher Mach simulations are under preparation.

5. High-Order Structure-Preserving Method Formulations for Compressible Turbulence in Nonuniform Curvilinear Grids

Selected formulations of high-order structure-preserving methods (SPM) in skew-symmetric splitting form for a curvilinear grid are discussed in this section. The discussion focuses on the Yee et al. and Sjogreen & Yee Entropy Split Methods that are entropy-conserving [19,23,26,30] since these entropy-conserving methods with recent development are not very familiar in the numerical method development community or among practitioners in the turbulence simulation community.

The entropy split method originally coined by Yee et al. in the late 1990s is the “Entropy Spitting Approach” [24,26].

The original entropy split method (the entropy-splitting approach) [26] is based on the Harten entropy function for the thermally perfect gas Euler equations. This original entropy splitting is a form of the skew-symmetric splitting of the nonlinear Euler flux derivatives [81] that takes advantage of the homogeneity property of the inviscid Euler flux and symmetrizable inviscid Euler flux derivatives. This splitting is of a form that is more suitable for the discrete stable energy norm estimate technique, including the summation-by-parts (SBP) boundary scheme estimate for an arbitrary order of non-dissipative central spatial schemes for non-periodic boundaries. It was later shown by the authors that this original entropy split method is entropy-conserving [18,19,21,22,29,76,94]. The Euler flux derivatives are approximated as a sum of a conservative portion and a non-conservative portion in conjunction with the SBP difference boundary closure [108] of Olsson & Oliger, Gerritsen & Olsson, and Yee et al. [26,88,92]. This original entropy split method is entropy-conserving, but the resulting method is not fully conservative. Two alternative methods were developed by the authors to overcome this partially conservative entropy split method. See [19,22,30] for details and test cases for a new class of entropy split methods that are entropy-conserving, and the resulting schemes are also conservative.

In 2022, instead of taking advantage of the homogeneity property of the inviscid Euler flux and symmetrizable inviscid Euler flux derivatives, a wider class of entropy split methods that do not require the homogeneous property was developed in [23,96]. The formulation is written in standard numerical flux form. This wider class of entropy split methods was extended to the equations of ideal MHD by the authors in [23,96].

For completeness, forms of other structure-preserving methods for comparison among the Yee et al. and Sjogreen & Yee entropy conserving, Tadmor-type of entropy conserving [82], momentum conserving [87], kinetic energy-preserving [86,89,90,91] methods and a combination of these structure-preserving methods [29,30,90,109] in a curvilinear grid are also included in this section.

The 3D compressible Euler equations of gas dynamics for a thermally perfect gas are a system of conservation laws,

$${\mathbf{q}}_t+{\mathbf{f}}_x+{\mathbf{g}}_y+{\mathbf{h}}_z=\mathbf{0}.$$

(26)

with conserved variables

$$\mathbf{q}={\left(\rho \rho u\rho v\rho we\right)}^T,$$

where $\rho$ denotes the density, $u,v,w$ are the velocities in the x-, y-, and z-directions, and e denotes the total energy.

When posed on a mapped domain with coordinate mapping $x=x(\xi ,\eta ,\zeta )$, $y=y(\xi ,\eta ,\zeta )$, $z=z(\xi ,\eta ,\zeta )$, the Equations (26) are transformed into

$$J{\mathbf{q}}_t+{(J{\xi }_x\mathbf{f}+J{\xi }_y\mathbf{g}+J{\xi }_z\mathbf{h})}_{\xi }+{(J{\eta }_x\mathbf{f}+J{\eta }_y\mathbf{g}+J{\eta }_z\mathbf{h})}_{\eta }+{(J{\zeta }_x\mathbf{f}+J{\zeta }_y\mathbf{g}+J{\zeta }_z\mathbf{h})}_{\zeta }=\mathbf{0}$$

on the unit cube, $(\xi ,\eta ,\zeta )\in {[0,1]}^3$, where J denotes the determinant of the Jacobian of the grid mapping.

The fluxes are considered in an arbitrary direction, $\mathbf{k}=\left(k_1{k}_2{k}_3\right)$, where, for mapped domains, $\mathbf{k}$ are the metric derivatives. For example, for the $\xi$-direction fluxes, $\mathbf{k}=\left(J{\xi }_{x}J{\xi }_{y}J{\xi }_z\right)$. The flux of the gas dynamics equation in the direction $\mathbf{k}$ is

$$\widehat{\mathbf{f}}=k_1\mathbf{f}+k_2\mathbf{g}+k_3\mathbf{h}={(\rho \widehat{u},\rho u\widehat{u}+k_{1}p,\rho v\widehat{u}+k_{2}p,\rho w\widehat{u}+k_{3}p,\widehat{u}(e+p))}^T.$$

(27)

Here, the velocity in the $\mathbf{k}$-direction is denoted as

$$\widehat{u}=k_{1}u+k_{2}v+k_{3}w.$$

The total energy is related to the pressure, p, via the ideal gas law,

$$e=\frac{p}{\gamma -1}+\frac{1}{2}\rho (u^2+v^2+w^2),$$

where $\gamma >1$ is a given constant. For simplicity, the 1D form of the equations is sometimes used in the description below. Writing x for the coordinate and $\mathbf{f}$ for the flux in one dimension should be interpreted for curvilinear grids as the coordinate $\xi$ with flux $\widehat{\mathbf{f}}$ and direction vector $\mathbf{k}$ equal to the corresponding metric derivatives.

Entropy is a convex function, $E=E\left(\mathbf{q}\right)$, such that smooth solutions of (26) satisfy the additional scalar conservation law

$$E_t+F_x+G_y+H_z=0,$$

(28)

where the entropy fluxes, F, G, and H, are related to the Euler fluxes via $E_{\mathbf{q}}{\mathbf{f}}_x=F_x$, and similarly for the y- and z-directions. The entropy variables are defined as

$$\mathbf{v}=E_{\mathbf{q}}\left(\mathbf{q}\right),$$

which is a well-defined change of variables, $\mathbf{v}=\mathbf{v}\left(\mathbf{q}\right)$, due to the convexity of $E\left(\mathbf{q}\right)$.

There are mainly two entropies for the Euler equations that are used in entropy-based numerical method developments: (a) Harten [81] considered the class of entropies

$$E_H=-\frac{\gamma +\alpha }{\gamma -1}\rho {\left(p{\rho }^{-\gamma }\right)}^{\frac{1}{\alpha +\gamma }},$$

(29)

where $\alpha$ is a parameter. To ensure that $E_H$ is convex, i.e., that the matrix ${\left(E_H\right)}_{\mathbf{q},\mathbf{q}}$ is positive-definite, $\alpha$ is required to satisfy $\alpha >0$ or $\alpha <-\gamma$. The corresponding x-direction entropy flux is $uE$. (b) The logarithmic entropy,

$$E_L=-\rho log\left(p{\rho }^{-\gamma }\right),$$

(30)

is another commonly used entropy. Entropy-conserving (EC) schemes are numerical discretizations for which the semi-discrete counterpart of (28) holds.

The high-order entropy split methods split the 3D thermally perfect gas dynamics Euler equations into conservative and non-conservative parts before discretization as

$${\mathbf{q}}_t+\frac{\beta }{\beta +1}({\mathbf{f}}_x^{\left(x\right)}+{\mathbf{f}}_y^{\left(y\right)}+{\mathbf{f}}_z^{\left(z\right)})+\frac{1}{\beta +1}(A^{\left(x\right)}{\mathbf{v}}_x+A^{\left(y\right)}{\mathbf{v}}_y+A^{\left(z\right)}{\mathbf{v}}_z)=\mathbf{0}.$$

(31)

The matrix in the x-direction is defined as

$$A^{\left(x\right)}={\mathbf{f}}_{\mathbf{v}}^{\left(x\right)}$$

(32)

with $\beta \ne -1$, and for physically relevant solutions, we must choose the positive split parameter $\beta >0$. The situation is similar for the other coordinate directions.

The original split form of the Euler flux derivative, but with negative $\beta$, was used in the original studies by Harten and Gerritsen & Olsson [81,88]. Yee et al., Vinokur & Yee, and Sjögreen et al. [24,26,36] found the physically relevant branch of the splitting parameter, $\beta >0$, by re-examining the results in [81,88]. Yee et al. [24,26] extended the entropy split form to include thermally perfect gas for moving curvilinear grids. In addition, they performed detailed numerical studies by applying high-order accurate entropy-splitting forms of the method and comparing the result with those of the standard unsplit form in various inviscid gas dynamics problems, including a 3D turbulent channel flow [16]. They found that using high-order central discretizations on the entropy split form of the Euler flux derivatives alone led to more stable computation than the standard unsplit approximation without additional numerical dissipation. Specifically, the long-time integration of selected smooth flows could be achieved with the entropy split scheme without the need for numerical dissipation, whereas numerical dissipation was required for the unsplit approximation. See [14,15,20,76,95,109,110] for some of the related entropy split methods’ developments. Comparisons of other forms of skew-symmetric splitting of the Euler flux derivatives [86,87,91] were also performed with results indicated in the aforementioned work.

Below, we expand the discussion of the original entropy split method using linear difference operators in one space dimension.

5.1. High-Order Entropy Split Methods Using Linear Difference Operators [19,20,22,26,76]

The original entropy split methods [26] that used linear difference operators for their construction were later proven to be entropy-conservative and stable for a thermally perfect gas in Sjögreen & Yee [19,20,22]. The various high-order methods resulting from applying classical spatial central, DRP, or Padé methods to the split form of the Euler flux derivative are entropy-conservative [18,26]. An entropy split approximation of the Euler flux derivatives based on DRP finite difference discretizations [75,93,111,112] was used in [18], and results for a aeroacoustic turbulence test case indicated a similar improvement in long-time integration numerical stability.

This original entropy split method construction differs from standard entropy-conserving methods in two ways. Firstly, the entropy split methods do not rely on a two-point numerical flux. Instead, the method is defined based on a standard linear difference operator. Secondly, the original entropy split method conserves entropy (29), but the method itself is not in a conservative form. See Sjögreen & Yee [22] for a recent new development to obtain entropy split methods that are entropy-conserving with the numerical method itself in conservation form.

The semi-discrete entropy split approximation of the one-dimensional version of (26) is

$$\frac{d}{dt}{\mathbf{q}}_j+\frac{\beta }{\beta +1}D{\mathbf{f}}_j+\frac{1}{\beta +1}{\left({\mathbf{f}}_{\mathbf{v}}\right)}_{j}D{\mathbf{v}}_j=0,j=1,\dots ,N,$$

(33)

where D is a linear finite difference operator (e.g., a high-order non-dissipative central difference, Padé, or DRP operator), and $\beta >0$ is a parameter related to $\alpha$ in (29) via $\beta =(\alpha +\gamma )/(1-\gamma )$.

The flux Jacobian matrix with respect to the entropy variables, ${\mathbf{f}}_{\mathbf{v}}$, is symmetric. In the most common case, D is a standard high-order SBP centered difference operator, but other operators are possible. For example, D could be a bandwidth-optimized operator with SBP closure, such as that developed in [75,93,111,112].

If the difference operator, D, has the SBP property, then the entropy-conserving property

$$\Delta x\frac{d}{dt}\sum _{j=1}^N{\omega }_{j}E\left(t\right)-F_1+F_N=0$$

(34)

holds; see [19].

Entropy splitting weights the non-conservative portion of the flux derivative via $1/(\beta +1)$, where $\beta >0$. Hence, using the range $\beta >0$ corresponds to giving the non-conservative portion a weight between zero and one, whereas $\beta <0$ leads, nonphysically, to a weight that is greater than one (for $-1<\beta <0$) or negative (for $\beta <-1$). Furthermore, because the non-conservative weight $1/(\beta +1)$ goes to zero as $\beta \to \infty$, choosing a large $\beta$ should reduce the influence of non-conservative effects.

In our early work performed on some smooth flow test cases using high-order non-dissipative discretizations, results indicated that $\beta =1$ (corresponding to a weight of $1/2$ on both the conservative portion and the non-conservative portion of the flux derivative) can improve stability without the need to add numerical dissipation. For shock-free turbulence, a similar improvement in stability was observed for $\beta =1$ or 2. For turbulence with shocklets, $\beta \ge 2$ is needed, possibly because larger $\beta$ gives better conservation properties. For turbulence with shocks, new forms of the entropy split method were proposed in Sjögreen & Yee [22] that include the derivation of a high-order conservative numerical flux for the non-conservative portion of the entropy splitting of the Euler flux derivatives in conjunction with their nonlinear filter approach.

5.2. Spatial Discretizations via Two-Point Numerical Fluxes

The majority of current high-order FDM methods are written in two-point numerical fluxes’ form. In this section, for simplicity of presentation, the numerical discretization of various methods is described in one space dimension using a uniform finite difference grid, $x_j=(j-1)\Delta x$, $j=1,\dots ,N$, where $\Delta x$ is the grid spacing.

The schemes will be considered on the semi-discrete form,

$$\frac{d{\mathbf{q}}_j}{dt}+\frac{{\mathbf{h}}_{j+1/2}-{\mathbf{h}}_{j-1/2}}{\Delta x}=\mathbf{0},$$

(35)

where ${\mathbf{h}}_{j+1/2}$ denotes numerical fluxes. For three space dimensions, the flux formulas are applied separately along each coordinate direction.

Here, we consider two-point second-order accurate fluxes, ${\mathbf{h}}_{j+1/2}=\mathbf{h}({\mathbf{q}}_j,{\mathbf{q}}_{j+1})$. These can be used as building blocks for higher-order methods, as is described in the next subsection.

The entropy-conserving property is [82],

$${({\mathbf{v}}_{j+1}-{\mathbf{v}}_j)}^T{\mathbf{h}}_{j+1/2}={\psi }_{j+1}-{\psi }_j,$$

(36)

where $\psi ={\mathbf{v}}^T\mathbf{f}-F$ is the entropy flux potential. Equation (36) is one equation for the five flux components, leading to non-unique EC methods. Interchanging j and $j+1$ in the differences on the right-hand-side and left-hand-side shows that $\mathbf{h}({\mathbf{q}}_{j+1},{\mathbf{q}}_j)$ also satisfies (36), making it reasonable to require the symmetry $\mathbf{h}({\mathbf{q}}_{j+1},{\mathbf{q}}_j)=\mathbf{h}({\mathbf{q}}_j,{\mathbf{q}}_{j+1})$. When (36) holds, it is possible to derive the semi-discrete entropy equality

$$\frac{dE\left({\mathbf{q}}_j\right)}{dt}+\frac{1}{\Delta x}(H_{j+1/2}-H_{j-1/2})=0$$

from (35), [82]. Here, $H_{j+1/2}$ are numerical entropy fluxes consistent with F.

The kinetic energy preservation (KEP) property of a three-point, second-order, accurate numerical semi-discretization holds if the semi-discrete solution satisfies

$$\frac{d{\left(e_k\right)}_j}{dt}+\frac{1}{\Delta x}(h_{j+1/2}^{\left(kin\right)}-h_{j-1/2}^{\left(kin\right)})+{\widehat{u}}_{j}D{p}_j=0$$

(37)

for any kinetic energy numerical fluxes, $h^{\left(kin\right)}$. The operator acting on the pressure is the standard second-order accurate $D{p}_j=(p_{j+1}-p_{j-1})/(2\Delta x)$. The simple identity $u_{j}D{p}_j=-p_{j}D{u}_j+(h_{j+1/2}^{\left(up\right)}-h_{j-1/2}^{\left(up\right)})$ with $h_{j+1/2}^{\left(up\right)}=(u_j{p}_{j+1}+u_{j+1}{p}_j)/2$ makes it possible to rewrite the non-conservative part of (37) as $p_{j}D{\widehat{u}}_j$.

The methods will be given in terms of their two-point, second-order, accurate numerical fluxes, ${\mathbf{h}}_{j+1/2}=\mathbf{h}({\mathbf{q}}_j,{\mathbf{q}}_{j+1})$.

From the two-point flux, a generalization to wider stencil methods with additional properties, such as higher-order accuracy or wavenumber optimization, is straightforward. Next follows a few results for generalized schemes.

5.3. Extension to High-Order Spatial Discretizations [20,76]

Definition 1.

Given a numerical flux function, ${\mathbf{h}}_{j+1/2}=\mathbf{h}({\mathbf{q}}_j,{\mathbf{q}}_{j+1})$, and a linear finite difference operator,

$${\mathbf{f}}_x\left(x_j\right)\approx D{f}_j=\frac{1}{\Delta x}\sum _{k=1}^p{\alpha }_k({\mathbf{f}}_{j+k}-{\mathbf{f}}_{j-k}),$$

(38)

the generalized numerical flux is

$${\mathbf{g}}_{j+1/2}=2\sum _{k=1}^p{\alpha }_k\sum _{m=1}^k\mathbf{h}({\mathbf{q}}_{j+m},{\mathbf{q}}_{j+m-k}).$$

(39)

Remark 4.

In addition, the generalized flux difference can be written:

$$\frac{1}{\Delta x}({\mathbf{g}}_{j+1/2}-{\mathbf{g}}_{j-1/2})=\frac{1}{\Delta x}\sum _{k=1}^{p}2{\alpha }_k(\mathbf{h}({\mathbf{q}}_{j+k},{\mathbf{q}}_j)-\mathbf{h}({\mathbf{q}}_j,{\mathbf{q}}_{j-k})).$$

(40)

Definition 2.

Given a numerical flux, ${\mathbf{h}}_{j+1/2}=\mathbf{h}({\mathbf{q}}_j,{\mathbf{q}}_{j+1})$, and a linear boundary operator of width q for the first s grid points,

$${\mathbf{f}}_x\left(x_j\right)\approx D_b{f}_j=\frac{1}{\Delta x}\sum _{k=1}^q{\beta }_{jk}{\mathbf{f}}_k,j=1,\dots ,s.$$

(41)

The generalized boundary operator is

$$\frac{1}{\Delta x}\sum _{k=1}^{q}2{\beta }_{jk}\mathbf{h}({\mathbf{q}}_k,{\mathbf{q}}_j),$$

(42)

Assuming symmetry, $\mathbf{h}({\mathbf{u}}_{j+1},{\mathbf{u}}_j)=\mathbf{h}({\mathbf{u}}_j,{\mathbf{u}}_{j+1})$, the following results hold.

Theorem 1.

If the linear operator (38) is pth-order accurate, then (40) is pth-order accurate.

Theorem 2.

If the two-point flux is EC, then (40) is EC. Furthermore, if the combined linear operator and boundary operator are SBP, then (42) and (40) satisfy the semi-discrete entropy balance

$$\Delta x\frac{d}{dt}\sum _{j=1}^N{\omega }_{j}E\left(t\right)-F_1+F_N=0,$$

where ${\omega }_j>0$ are SBP norm weights, and $F_1,F_N$ are boundary entropy fluxes.

If the linear operator is wavenumber-optimized (e.g., DRP), then (40) is a generalization of the wavenumber-optimized operator. However, since the wavenumber optimization is a linear concept, it is harder to define a wavenumber property of the generalized flux, other than noting that it becomes the wavenumber-optimized method wherever the flux is linear. However, Theorem 2 shows that the generalized wavenumber-optimized approximation is entropy-conserving.

For the generalized scheme, the KEP property is defined using (37), but the operator, D, acting on the pressure now is required to be (38).

5.4. High-Order Numerical Fluxes in Curvilinear Grids for Comparison Study [20,76]

Denote the arithmetic average of a grid function, u, over two grid points using

$$\left\{u\right\}=(u_{j+1}+u_j)/2,$$

where $u_j$ denotes the numerical approximation of variable u at grid point $x_j$.

The five flux components of the Euler equations will be denoted using the superscript $\left(m\right),m=1,\dots ,5$; e.g., the mass flux is denoted as $h^{\left(1\right)}$.

5.4.1. Kinetic Energy-Preserving (KEP) Property

A simple way to verify a numerical method through KEP is the following.

Theorem 3.

The KEP holds if the numerical momentum fluxes satisfy

$$\begin{array}{cc}\hfill h^{\left(2\right)}& =\left\{u\right\}{h}^{\left(1\right)}+k_1\left\{p\right\}\hfill \end{array}$$

(43)

$$\begin{array}{cc}\hfill h^{\left(3\right)}& =\left\{v\right\}{h}^{\left(1\right)}+k_2\left\{p\right\}\hfill \end{array}$$

(44)

$$\begin{array}{cc}\hfill h^{\left(4\right)}& =\left\{w\right\}{h}^{\left(1\right)}+k_3\left\{p\right\}.\hfill \end{array}$$

(45)

See [90].

Proof. 

The identity

$$\frac{1}{2}{\left(\rho u^2\right)}_t=u{\left(\rho u\right)}_t-\frac{1}{2}{u}^2{\left(\rho \right)}_t$$

shows that

$$\frac{d{\left(e_k\right)}_j}{dt}=u_j\frac{d{\rho }_j{u}_j}{dt}+v_j\frac{d{\rho }_j{v}_j}{dt}+w_j\frac{d{\rho }_j{w}_j}{dt}-\frac{1}{2}(u_j^2+v_j^2+w_j^2)\frac{d{\rho }_j}{dt}.$$

Using (35) to replace time derivatives with spatial differences gives

$$\begin{array}{c}\Delta x\frac{d{\left(e_k\right)}_j}{dt}=u_j(h_{j+1/2}^{\left(2\right)}-h_{j-1/2}^{\left(2\right)})-\frac{1}{2}{u}_j^2(h_{j+1/2}^{\left(1\right)}-h_{j-1/2}^{\left(1\right)})\hfill \\ +v_j(h_{j+1/2}^{\left(3\right)}-h_{j-1/2}^{\left(3\right)})-\frac{1}{2}{v}_j^2(h_{j+1/2}^{\left(1\right)}-h_{j-1/2}^{\left(1\right)})\\ \hfill +w_j(h_{j+1/2}^{\left(4\right)}-h_{j-1/2}^{\left(4\right)})-\frac{1}{2}{w}_j^2(h_{j+1/2}^{\left(1\right)}-h_{j-1/2}^{\left(1\right)}).\end{array}$$

(46)

If (43)–(45) hold, then

$$\begin{array}{c}{u}_j(h_{j+1/2}^{\left(2\right)}-h_{j-1/2}^{\left(2\right)})-\frac{1}{2}{u}_j^2(h_{j+1/2}^{\left(1\right)}-h_{j-1/2}^{\left(1\right)})\hfill \\ =\frac{1}{2}(u_j(u_{j+1}+u_j)-u_j^2)h_{j+1/2}^{\left(1\right)}-\frac{1}{2}(u_j(u_j+u_{j-1})-u_j^2)h_{j-1/2}^{\left(1\right)}+k_1{u}_j\frac{1}{2}(p_{j+1}+p_j-(p_j+p_{j-1}))\hfill \\ =\frac{1}{2}(u_{j+1}{u}_j{h}_{j+1/2}^{\left(1\right)}-u_j{u}_{j-1}{h}_{j-1/2}^{\left(1\right)})+k_1{u}_j\frac{1}{2}(p_{j+1}-p_{j-1}).\hfill \end{array}$$

(47)

The terms multiplied by $v_j$ and $w_j$ in (46) are rewritten in the same way to obtain (37) with

$$h_{j+1/2}^{\left(kin\right)}=\frac{1}{2}(u_{j+1}{u}_j+v_{j+1}{v}_j+w_{j+1}{w}_j)h_{j+1/2}^{\left(1\right)}.$$

If the KEP property is hard to verify analytically, it is possible to apply the scheme to a flow that satisfies periodic boundary conditions and check numerically whether

$$\sum _{j=1}^{N-1}{u}_j{r}_j^{\left(2\right)}+v_j{r}_j^{\left(3\right)}+w_j{r}_j^{\left(4\right)}-\frac{1}{2}(u_j^2+v_j^2+w_j^2)r_j^{\left(1\right)}=\sum _{j=1}^{N-1}\widehat{u_j}D{p}_j$$

(48)

holds. Here, $r_j^{\left(m\right)}$ denotes the flux difference $(h_{j+1/2}^{\left(m\right)}-h_{j-1/2}^{\left(m\right)})/\Delta x$. If a numerical example can be found where (48) does not hold, then clearly, KEP cannot hold. If (48) seems to hold numerically, the test is inconclusive.

Below, a description of a few methods, given in terms of their two-point numerical flux functions, is provided. □

5.4.2. Tadmor-Type Entropy-Conserving Numerical Flux Using the Entropy Function $E_L=-\rho log\left(p{\rho }^{-\gamma }\right)$: ECLOG

The Tadmor-type entropy-conserving method conserves entropy (30). In its original form, it is not a kinetic energy preservation. The logarithmic average of the density is defined via

$${\rho }^{ln}=\frac{{\rho }_{j+1}-{\rho }_j}{log{\rho }_{j+1}-log{\rho }_j}.$$

(49)

The numerical flux vector is

$$\begin{array}{cc}\hfill h^{\left(1\right)}& ={\rho }^{ln}\left\{\widehat{u}\right\}\hfill \end{array}$$

(50)

$$\begin{array}{cc}\hfill h^{\left(2\right)}& ={\rho }^{ln}\left\{\widehat{u}\right\}\left\{u\right\}+k_1\frac{\left\{\rho \right\}}{\{\rho /p\}}\hfill \end{array}$$

(51)

$$\begin{array}{cc}\hfill h^{\left(3\right)}& ={\rho }^{ln}\left\{\widehat{u}\right\}\left\{v\right\}+k_2\frac{\left\{\rho \right\}}{\{\rho /p\}}\hfill \end{array}$$

(52)

$$\begin{array}{cc}\hfill h^{\left(4\right)}& ={\rho }^{ln}\left\{\widehat{u}\right\}\left\{w\right\}+k_3\frac{\left\{\rho \right\}}{\{\rho /p\}}\hfill \end{array}$$

(53)

and

$$\begin{array}{c}{h}^{\left(5\right)}=\frac{1}{\gamma -1}\frac{{\rho }^{ln}}{{(\rho /p)}^{ln}}\left\{\widehat{u}\right\}+\frac{\left\{\rho \right\}}{\{\rho /p\}}\left\{\widehat{u}\right\}-\frac{1}{2}{\rho }^{ln}\left\{\widehat{u}\right\}(\left\{u^2\right\}+\left\{v^2\right\}+\left\{w^2\right\})\hfill \\ \hfill +{\rho }^{ln}\left\{\widehat{u}\right\}({\left\{u\right\}}^2+{\left\{v\right\}}^2+{\left\{w\right\}}^2).\end{array}$$

(54)

Then, the logarithmic average of $\rho /p$ required in the energy flux is defined via (49) with $\rho /p$ replacing $\rho$.

5.4.3. Entropy-Conserving Numerical Flux of Tadmor-Type Using the Entropy Function $E_L=-\rho log\left(p{\rho }^{-\gamma }\right)$ with the KEP Property [90]: ECLOGKP

This Tadmor-type method conserves entropy (30). Using the Ranocha formulation, it has the kinetic energy preservation property.

For this entropy-conserving and kinetic energy-preserving method, the numerical flux vector is defined via

$$\begin{array}{cc}\hfill h^{\left(1\right)}& ={\rho }^{ln}\left\{\widehat{u}\right\}\hfill \end{array}$$

(55)

$$\begin{array}{cc}\hfill h^{\left(2\right)}& ={\rho }^{ln}\left\{\widehat{u}\right\}\left\{u\right\}+k_1\left\{p\right\}\hfill \end{array}$$

(56)

$$\begin{array}{cc}\hfill h^{\left(3\right)}& ={\rho }^{ln}\left\{\widehat{u}\right\}\left\{v\right\}+k_2\left\{p\right\}\hfill \end{array}$$

(57)

$$\begin{array}{cc}\hfill h^{\left(4\right)}& ={\rho }^{ln}\left\{\widehat{u}\right\}\left\{w\right\}+k_3\left\{p\right\}\hfill \end{array}$$

(58)

and

$$\begin{array}{c}{h}^{\left(5\right)}=\frac{1}{\gamma -1}\frac{{\rho }^{ln}}{{(\rho /p)}^{ln}}\left\{\widehat{u}\right\}+\left\{p\right\}\left\{\widehat{u}\right\}-\frac{1}{2}{\rho }^{ln}\left\{\widehat{u}\right\}(\left\{u^2\right\}+\left\{v^2\right\}+\left\{w^2\right\})\hfill \\ \hfill +{\rho }^{ln}\left\{\widehat{u}\right\}({\left\{u\right\}}^2+{\left\{v\right\}}^2+{\left\{w\right\}}^2)-\frac{1}{4}(p_{j+1}-p_j)({\widehat{u}}_{j+1}-{\widehat{u}}_j).\end{array}$$

(59)

5.4.4. Tadmor-Type Entropy-Conserving Numerical Flux Using the Harten Entropy Function $E_H=-\frac{\gamma +\alpha }{\gamma -1}\rho {\left(p{\rho }^{-\gamma }\right)}^{\frac{1}{\alpha +\gamma }}$ with the KEP Property [90]: ECHKP

This method conserves entropy (29), as well as automatically possessing the KEP property. Denote

$$z=\frac{\rho }{p}{\left(p{\rho }^{-\gamma }\right)}^{\frac{1}{\alpha +\gamma }}.$$

(60)

The numerical two-point flux is defined as

$$\begin{array}{cc}\hfill h^{\left(1\right)}& =\frac{\left\{z\widehat{u}\right\}}{\left\{z^{-\gamma /\alpha }\right\}{\left\{p^{(1-\gamma -\alpha )/\alpha }\right\}}_e}\hfill \\ \hfill h^{\left(2\right)}& =h^{\left(1\right)}\left\{u\right\}+k_1\left\{p\right\}\hfill \\ \hfill h^{\left(3\right)}& =h^{\left(1\right)}\left\{v\right\}+k_2\left\{p\right\}\hfill \\ \hfill h^{\left(4\right)}& =h^{\left(1\right)}\left\{w\right\}+k_3\left\{p\right\}\hfill \end{array}$$

$$\begin{array}{c}\hfill h^{\left(5\right)}=h^{\left(1\right)}(\frac{\gamma }{\gamma -1}\left\{p^{-(\gamma -1)/\alpha }\right\}{\left\{z^{-\frac{\gamma +\alpha }{\alpha }}\right\}}_e+{\left\{u\right\}}^2+{\left\{v\right\}}^2+{\left\{w\right\}}^2-\frac{1}{2}\{u^2+v^2+w^2\}).\end{array}$$

The exponential average is defined via

$${\left\{q^{\beta }\right\}}_e=(q_{j+1}^{\beta +1}-q_j^{\beta +1})/\left((\beta +1)(q_{j+1}-q_j)\right).$$

(61)

5.4.5. Momentum-Conserving Method of Ducros et al. [19,87]: DS

For the Ducros et al. momentum-conserving method, the corresponding numerical flux vector is defined as

$$\begin{array}{cc}\hfill h^{\left(1\right)}& =\left\{\rho \right\}\left\{\widehat{u}\right\}\hfill \end{array}$$

(62)

$$\begin{array}{cc}\hfill h^{\left(2\right)}& =\left\{\rho u\right\}\left\{\widehat{u}\right\}+k_1\left\{p\right\}\hfill \end{array}$$

(63)

$$\begin{array}{cc}\hfill h^{\left(3\right)}& =\left\{\rho v\right\}\left\{\widehat{u}\right\}+k_2\left\{p\right\}\hfill \end{array}$$

(64)

$$\begin{array}{cc}\hfill h^{\left(4\right)}& =\left\{\rho w\right\}\left\{\widehat{u}\right\}+k_3\left\{p\right\}\hfill \end{array}$$

(65)

$$\begin{array}{cc}\hfill h^{\left(5\right)}& =\{e+p\}\left\{\widehat{u}\right\}.\hfill \end{array}$$

(66)

Remark 5.

The identity

$$\frac{1}{2}D\left(ab\right)+\frac{1}{2}aD\left(b\right)+\frac{1}{2}bD\left(a\right)=\frac{1}{4\Delta x}\left((a_{j+1}+a_j)(b_{j+1}+b_j)-(a_j+a_{j-1})(b_j+b_{j-1})\right)$$

with $Da=(a_{j+1}-a_{j-1})/2\Delta x$ shows that the two-factor fluxes can be interpreted as a conservative/non-conservative split approximation. By rewriting the fluxes, they can be written in a conservation form. See [19] for the generalization of the Ducros et al. momentum-conserving method to an arbitrary order of accuracy in a conservative form. From our numerical comparison, this method turned out to be fairly robust in computations, but it is not known to conserve any entropy.

5.4.6. Ducros et al. Momentum-Conserving Method with KEP Property [90]: DSKP

The DSKP method is a variant of DS that can be easily formulated to include the Ranocha KEP property.

The corresponding numerical flux vector in this case is defined as

$$\begin{array}{cc}\hfill h^{\left(1\right)}& =\left\{\rho \widehat{u}\right\}\hfill \end{array}$$

(67)

$$\begin{array}{cc}\hfill h^{\left(2\right)}& =\left\{\rho \widehat{u}\right\}\left\{u\right\}+k_1\left\{p\right\}\hfill \end{array}$$

(68)

$$\begin{array}{cc}\hfill h^{\left(3\right)}& =\left\{\rho \widehat{u}\right\}\left\{v\right\}+k_2\left\{p\right\}\hfill \end{array}$$

(69)

$$\begin{array}{cc}\hfill h^{\left(4\right)}& =\left\{\rho \widehat{u}\right\}\left\{w\right\}+k_3\left\{p\right\}\hfill \end{array}$$

(70)

$$\begin{array}{cc}\hfill h^{\left(5\right)}& =\left\{\rho \widehat{u}\right\}\left\{\frac{e+p}{\rho }\right\}.\hfill \end{array}$$

(71)

5.4.7. Kennedy–Gruber–Pirozzoli Kinetic Energy-Preserving Method: KGP

The Kennedy–Gruber–Pirozzoli kinetic energy-preserving method [86,91] turned out to perform very well among the considered methods and test case comparisons. The numerical flux vector is defined as

$$\begin{array}{cc}\hfill h^{\left(1\right)}& =\left\{\rho \right\}\left\{\widehat{u}\right\}\hfill \end{array}$$

(72)

$$\begin{array}{cc}\hfill h^{\left(2\right)}& =\left\{\rho \right\}\left\{\widehat{u}\right\}\left\{u\right\}+k_1\left\{p\right\}\hfill \end{array}$$

(73)

$$\begin{array}{cc}\hfill h^{\left(3\right)}& =\left\{\rho \right\}\left\{\widehat{u}\right\}\left\{v\right\}+k_2\left\{p\right\}\hfill \end{array}$$

(74)

$$\begin{array}{cc}\hfill h^{\left(4\right)}& =\left\{\rho \right\}\left\{\widehat{u}\right\}\left\{w\right\}+k_3\left\{p\right\}\hfill \end{array}$$

(75)

$$\begin{array}{cc}\hfill h^{\left(5\right)}& =\left\{\rho \right\}\left\{\widehat{u}\right\}\{(e+p)/\rho \}.\hfill \end{array}$$

(76)

Remark 6.

Similarly to the Ducros splitting, KGP splitting can be written as a conservative/non-conservative split approximation, but using three factors instead of two:

$$\begin{array}{c}\frac{1}{4}(D\left(abc\right)+abD\left(c\right)+acD\left(b\right)+bcD\left(a\right)+aD\left(bc\right)+bD\left(ac\right)+cD\left(ab\right))\hfill \\ \hfill =\frac{1}{8\Delta x}\left((a_{j+1}+a_j)(b_{j+1}+b_j)(c_{j+1}+c_j)-(a_j+a_{j-1})(b_j+b_{j-1})(c_j+c_{j-1})\right).\end{array}$$

(77)

6. Entropy Split Method and Recent Extensions ES

For the extension of the original entropy split method, the semi-discrete entropy split approximation of the 1D version of (26),

$$\frac{d}{dt}{\mathbf{q}}_j+\frac{\beta }{\beta +1}D{\mathbf{f}}_j+\frac{1}{\beta +1}{\left({\mathbf{f}}_{\mathbf{v}}\right)}_{j}D{\mathbf{v}}_j=0,j=1,\dots ,N,$$

(78)

is repeated here, where D is a linear finite difference operator (38), and $\beta >0$ is a parameter related to $\alpha$ in (29) via $\beta =(\alpha +\gamma )/(1-\gamma )$.

The flux Jacobian matrix with respect to the entropy variables, ${\mathbf{f}}_{\mathbf{v}}$, is symmetric. In the most common case, D is a standard high-order SBP-centered difference operator, but other operators are possible. For example, D could be a bandwidth-optimized operator with SBP closure, such as that developed in [93].

If the difference operator, D, has the SBP property, then the entropy-conserving property

$$\Delta x\frac{d}{dt}\sum _{j=1}^N{\omega }_{j}E\left(t\right)-F_1+F_N=0$$

can be proven for (78); see [19].

It is noted that, for $\beta =1$, the entropy split method using a central spatial scheme, after being rewritten into conservative differencing, becomes the standard Ducros et al. splitting, and it can easily be rewritten into a conservation form.

6.1. Hybridized Entropy Split Method with Ducros et al. Splitting: ESDS

For this method, we hybridized DS with ES by replacing the centered approximation of the conservative part of ES with the DS flux difference:

$$\frac{d}{dt}{\mathbf{q}}_j+\frac{\beta }{\beta +1}\frac{1}{\Delta x}({\mathbf{h}}_{j+1/2}-{\mathbf{h}}_{j-1/2})+\frac{1}{\beta +1}{\left({\mathbf{f}}_{\mathbf{v}}\right)}_{j}D{\mathbf{v}}_j=0,j=1,\dots ,N,$$

(79)

where ${\mathbf{h}}_{j+1/2}$ is the DS numerical flux function. The performance of the ESDS method is shown in the comparison of high-order methods in the next section.

6.2. New Conservative Entropy Split Method ESSW

To overcome the issue of being a non-conservative discretization via the original entropy-conserving split scheme, Sjögreen & Yee [22] derived a high-order conservative numerical flux for the non-conservative portion of the entropy splitting of the Euler flux derivatives. Due to the construction, this conservative numerical flux requires a higher operations count, and it is less stable than the original semi-conservative split method. For the derivation, formulation, and test case comparison of the scheme NEW, see Sjögreen & Yee [22]. The comparison indicated that the Tadmor entropy-conservative method [82] of the same order requires a higher operations count than the NEW entropy split construction. Aside from being entropy-stable, the main advantage of the entropy split method is its superior stability without added numerical dissipation over the standard centered difference method. In addition, the entropy split method is more stable for the long-time integration of smooth flows than the conservative Ducros et al. splitting [87] or the Kennedy–Gruber–Pirozzoli [86,89] splitting. Its drawback is that it is non-conservative. When the entropy split method is used (together with conservative shock-capturing dissipation as a nonlinear filter step) to compute flows with shock waves, the degree of incorrect shock speed increases as the shock strength increases.

An alternative to using the conservative version of the entropy split form for problems containing shock waves is to use a shock detector to switch the entropy split discretization to a standard centered approximation in the neighborhood of shock waves to avoid the wrong shock speed. Alternatively, since entropy splitting is more stable for the long-time integration of smooth flows than the conservative Ducros et al. splitting [87] or the Kennedy–Gruber–Pirozzoli [86,89] splitting, but the latter two splittings are more stable than the pure central spatial discretization, the switch can be between the entropy split method and conservative skew-symmetric splitting, e.g., Ducros et al. or Kennedy–Gruber–Pirozzoli splitting.

An extension of these structure-preserving methods to the ideal MHD equations can be found in our previous work [20,30]. For the extension of the new entropy split formulation using a two-point numerical flux without requiring the homogeneity property of the compressible inviscid flux, see the last section or [23]. Since this new approach was just published, a large portion of [23] is repeated in Section 8.

7. Gas Dynamics Comparison Among High-Order Methods

This section illustrates different test cases compared with selected high-order structure-preserving methods before a discussion of the extension of the new entropy split method for the equation of ideal MHD. The same notation and nomenclature for these methods indicated in our published work are defined below.

7.1. Spatial Discretization with Various Structure-Preserving Properties

The nomenclature for the high-order structure-preserving methods considered for test cases’ comparison is as follows:

• ECHKP: Entropy conserving using the Harten class of entropy functions, $E_H=-\frac{\gamma +\alpha }{\gamma -1}\rho {\left(p{\rho }^{-\gamma }\right)}^{\frac{1}{\alpha +\gamma }}$ [81,82]. It turned out that this method, in its base form, also satisfies Ranocha’s kinetic energy preservation condition (KEP); so, there is only one variant of this method [76,90].
• ECLOG: The Tadmor-type entropy-conserving method using the Tadmor entropy function $E_L=-\rho log\left(p{\rho }^{-\gamma }\right)$.
• ES: The skew-symmetric splitting of the inviscid flux derivative that is entropy-conserving and entropy-stable using the Harten entropy function [81] and the generalized energy norm with summation by parts (SBP) [19,26,94].
• DS: A momentum-conserving Ducros et al. skew-symmetric split of the inviscid flux derivative [87].
• KGP: Kennedy–Gruber–Pirozzoli (KGP) skew-symmetric splitting of the inviscid flux derivative that is kinetic energy-preserving [86,89,91].
• ESDS: Entropy split with Ducros et al. splitting [19].
• ESSW: Entropy split with Ducros et al. splitting but a switch to regular central near-discontinuities [22].
• ECLOGKP: A Tadmor-type entropy-conserving method using the Tadmor entropy function with Ranocha’s kinetic energy-preserving modification [90].
• DSKP: Ducros split with KEP.

According to our previous studies [19,20,22], the Tadmor-type entropy conserving methods EC, ECLOG, and ECLOGKP are the most CPU-intensive methods among the nine methods. For the results by EC, see [19,20,22]. It uses approximately twice the CPU per time step as the ES, ESDS, and ESSW methods. DS is the least CPU-intensive. A comparison of execution times was given in previous published works.

For more test case comparisons, including strong shock waves, see [19,20,22]. Here, only the 2D isotropic vortex convection, 3D Taylor–Green vortex shock-free turbulence, and the Brio–Wu MHD shock tube test cases were selected for the nine-method comparison. The comparison included the maximum norm error, mass conservation errors, entropy, and kinetic energy vs. the time for the eighth-order methods. It is noted that the performance comparison discussed here pertained to the chosen flow type, governing equation set, and uniform grid spacings without any grid adaptation. The performance of the entropy split methods as a function of the split parameter $\beta$ for other test cases can be found in our previously published work indicated above. Just like other current high-order method developments in the literature, the performance of the entropy split method is highly dependent on the grid size, flow type, flow condition, shock-free turbulence, turbulence with shocklets, and turbulence with strong shocks. Below, only results by the eighth-order classical spatial discretization are shown. Studies using DRP and compact spatial discretizations were also performed, but they are not shown here. See [19,20,21,22,29] for additional investigations.

7.2. Smooth Flow Test Case for Gas Dynamics: 2D Isotropic Vortex Convection the Popular Test Case to Test the Stability and Accuracy of Long-Time Integration for a Smooth Flow Is the 2D Isentropic Vortex Convection with the Initial Data Indicated in Figure 24

The exact solution is the initial data translated, $\mathbf{u}(x,t)={\mathbf{u}}_0(x-u_{\infty }t,y-v_{\infty }t)$. The computational domain was of the size $0\le x\le 18$, $0\le y\le 18$ with periodic boundary conditions. The center of the vortex was $(x_0,y_0)=(9,9)$. We integrated this test case to a final time $T=800$ for a majority of the comparisons. Two grids, ${101}^2$ and ${201}^2$, were considered for the comparison. For $\beta =2$ and the coarse grid ${101}^2$, we integrated a final time of $T=1440$, which is more than 100 times longer than that reported in the literature and twice as long as that in most of our previous studies with an end time of 720. To obtain better time accuracy, the CFL number was set to $0.4$, leading to a fairly small time step.

Investigating the level of the maximum norm error via the eighth-order ES and ESSW methods as a function of $\beta$ vs. time using the two grids, Figure 25 shows that a lower maximum norm error can be prolonged for longer time integration using $\beta =1$ than $\beta =2$ with the finer grid for both the ES and ESSW methods. For the same grid, the maximum norm error is highly dependent on the value of $\beta$. The error remains the same for $\beta =1$ in both the ES and ESSW methods. However, for other $\beta =2$ and $2.5$, the error using the ES method is slightly different from that of the ESSW method.

Figure 24. 2D isentropic vortex convection problem setup.

Figure 24. 2D isentropic vortex convection problem setup.

Figure 26 shows the maximum norm error, mass conservation errors, entropy, and kinetic energy for the nine eighth-order methods using $\beta =1$ on the fine grid. Overall, ES performed the best except in the conservation of mass breaks down at around time 680. The kinetic energy and entropy results show the quantity with its value at time zero subtracted; e.g., the kinetic energy ($E_{kin\left(t\right)}$) shown is $E_{kin}\left(t\right)-E_{kin}\left(0\right)$.

To show the performance of $\beta =2$, Figure 27 shows the time evolution to the final end time of $T=1440$ using the coarse grid ${100}^2$ in terms of the maximum norm error, the two entropies, and the kinetic energy as a function of time for the eight methods using $\beta =2$ in the entropy split methods. The performance of the ES method using $\beta =2$ was very different from that using $\beta =1$ with the same grid.

The entropy $E_L$ shows the quantity with its value at time zero subtracted. The KGS, DS, and ESDS methods blew up before time 300, while the four methods that conserve entropy all ran to the final time. However, the methods that did not blow up had very large errors after time 300.

The entropies were not perfectly conserved for longer times. The likely explanation is that poor accuracy for long times made the time discretization error significant. This, in turn, destroyed the conservation when entropy conservation held only for the semi-discrete problem. However, as seen from the zoomed-in plots, the entropy and kinetic energy were very close to being constant up to the time when the errors became large. This is an improvement that we expected to see with the more accurate time integration.

The results for ES and ESDS highly depended on the split parameter $\beta$. Only comparisons using $\beta =1$ and $\beta =2$ based on the study in [19] are shown here. Using $\beta =1$ with grid refinement, the performance of ES was stable and more accurate for a longer time integration. In summary, for the 2D vortex convection case, the ESDS, ESSW, and KGP methods conserve entropy over a shorter time integration and are less accurate for a longer time integration than the other considered methods. For ESDS and ESSW, their modifications from the original ES method might interfere with the entropy conservation property.

For the splitting parameter $\beta$, even for shock-free long-time integration flows, for physical argument, it is preferred to use $\beta =1$ or greater, as, for $\beta <1$, more than a 50% non-conservative portion of the split flux derivative is used.

7.3. Three-Dimensional Shock-Free Compressible Turbulence Gas Dynamics Test Case—3D Taylor–Green Vortex

The well-known shock-free compressible turbulence test case to evaluate the stability and accuracy for gas dynamics is the Taylor–Green vortex [113]. The 3D Euler equations of compressible gas dynamics are solved with $\gamma =5/3$. The computational domain is a cube with sides of length $2\pi$ and with periodic boundary conditions in all three directions. The initial conditions are

$$\begin{array}{cc}\hfill & \rho =1p=100+\left\{[cos(2z)+2)(cos(2x)+cos(2y\left)\right]-2\right\}/16\hfill \end{array}$$

(80)

$$\begin{array}{cc}\hfill & u=sinxcosycosz,v=-cosxsinycosz,w=0.\hfill \end{array}$$

(81)

The problem is solved to time 20. In our previous studies, the uniform coarse grid DNS of ${64}^3$ grid points solutions was compared with the filtered fine grid DNS solutions using a uniform ${256}^3$ grid points. Here, we used the same uniform coarse grid to examine the nonlinear stability and accuracy of the three new eight-order NEW, ESSW, and ESDS methods. The total kinetic energy of the exact solution is constant in time.

The present discussion of numerical results was confined to a coarse-grid DNS comparison among the methods. It is noted that, for this Taylor–Green inviscid problem, small scales were generated that eventually caused large errors in the solution due to an inadequate resolution. This probably occurred around $T=5$. Another issue is that, for very low dissipative or non-dissipative numerical methods for the simulation of turbulent flows, even with extreme grid refinement, grid convergence cannot be obtained, as the original inviscid Euler equations are chaotic in nature. With sufficient but not excessive numerical dissipations, one is solving the equivalent of a very high Reynolds number in Navier–Stokes equations. See Yee & Sjögreen [35] for a study. The end time is 20 instead of the standard end time, 10, so that the solution behavior could be observed twice when using the same RK4 time discretization and CFL number, 0.4. Harten’s entropy was used for all schemes except ECLOG and ECLOGKP, for which we used the log entropy. Furthermore, the entropy $E\left(t\right)$ was normalized to show the deviation from the initial entropy. Thus, $E\left(t\right)-E\left(0\right)$ is shown in the figures.

Figure 28 shows the comparison of kinetic energy, the maximum norm error, and entropy vs. the time for the nine eighth-order methods. It is interesting to see the behavior of doubling the time integration duration for such coarse-grid DNS computations.

The ESDS method became unstable at around time 6, and the DSKP method became unstable at around time 7. All other schemes ran to completion. Again, the kinetic energy and entropy results show the quantity with its value at time zero subtracted; e.g., the kinetic energy ($E_{kin\left(t\right)}$) shown is $E_{kin}\left(t\right)-E_{kin}\left(0\right)$. The ES method started to lose some energy at a later time. Otherwise, the stable results were similar. ECLOGKP and KGS are indistinguishable, and they fall on top of each other in the zoomed-in figure. One surprising result is that the ECBKP method was expected to preserve the kinetic energy in the same manner but did not do so. The other stable methods were a little off, but this is only visible in the close-up. The ES and ECBKP methods fell on top of each other, as we would expect, since these two schemes conserve entropy in a discretized sense. The ECLOGKP method was also on top of ES and ECBKP, making it hard to visualize the differences in the results.

Not included here are studies illustrating that the logarithmic entropy function conserved Harten’s entropy almost perfectly. The ECLOGKP and ECLOG methods conserve entropy, as illustrated in the figure where their solutions are shown on top of each other. Overall, ECLOGKP, ECLOG, ECBKP, KGS, and DS are very similar. One has to zoom in very much in this region to see any differences for this test case. However, differences might be larger for other flow problems. As can be seen in this test case, the ES method behaves somewhat differently. DSKP and ESDS do not perform well. It is noted again that the results for ES and ESDS are highly dependent on $\beta$. Only results with $\beta =2$ are shown, based on the study in [19].

Figure 29 compares how well these eighth-order methods conserve mass and momentum in the x-, y-, and z-directions. The figures show the total mass (or momentum) as a function of time, normalized so that the initial mass (or momentum) is zero. For the mass-conservation plot, all schemes except ES, ESSW, and ESDS are perfectly conservative. The normalized $M\left(t\right)-M\left(0\right)$ is zero as long as the method is stable for all methods except the ES scheme; these values are around $e^{-9}$ most of the time. In addition, these figures are shown in the log scale, for which we artificially set the schemes with zero conservation error to $e^{-16}$.

The Taylor–Green vortex test case has a similar method comparison conclusion as the vortex convection test case. The ESDS, ESSW, and KGP methods conserve entropy in a shorter time integration and are less accurate for a longer time integration than the other considered methods. For ESDS and ESSW, their modification from the original ES method might interfere with the entropy conservation property.

7.4. A Difficult Brio–Wu MHD Shock Tube Test Case

For the extension of these high-order methods for the ideal MHD in curvilinear grids, see [30].

Here, we include the performance of a very difficult Brio–Wu MHD shock tube problem. It involves two fast rarefaction waves, a slow compound wave, contact discontinuity, and a slow shock wave. According to the Flash code [45,114], standard high-order shock-capturing methods exhibit oscillatory solutions. Examples are PPM (the piecewise parabolic method of Collela & Woodward [115]) and third- to fifth-order WENO (WENO3–WENO5). Most researchers resorted to using first-order or very diffusive second-order methods with local Lax–Freidrichs (LLF) or variants of the HLL-type numerical fluxes [116,117]. The oscillations increase with an increase in $B_y$; the reason is that a stronger $B_y$ introduces more of a transverse effect that resists shock propagation in the x-direction, causing the shock to move slowly.

The test case was solved on domain $0\le x\le 1$ to time $0.1$, using 801 uniformly distributed grid points. A grid refinement study was also performed using a 1601 uniform grid.

Details of the test case can be found in [23,30]. Figure 30 shows the problem setup. Figure 31 shows the comparison among seven eighth-order structure-preserving base schemes using the same WENO7 as the nonlinear filter with the flow sensor indicated in [30]. The figure shows the comparison of density and pressure. From the figure, one can observe the different resolutions at the discontinuity locations. The results demonstrate that all the considered high-order structure-preserving methods perform well on a rapidly developing flow consisting of shocks and rarefaction waves but without complex flow features. The use of structure preservation high-order methods in this case only helps with numerical stability but offers less of an improvement in accuracy. However, these high-order nonlinear filter structure-preserving methods show a great improvement compared to the results indicated in [45,114,116,117].

7.5. Summary

In summary, for the chosen test cases, the long-time integration performance of the entropy split methods as a function of the split parameter $\beta$ is dependent on the flow type, flow condition, and grid spacing. For finer grids, the stability of long-time integration for the same $\beta$ value can be different. The general guideline is that, for very smooth pure convection flows, $\beta =1$ (50% conservative and 50% non-conservative of the split term) or a $\beta$ between 1 and 2 can improve the stability tremendously without added numerical dissipation. A finer grid in general prolongs stability for a longer time integration. The split method conserves mass and entropy. For problems with shock waves, a $\beta$ between 20 and 26 performs well with the NEW and ESSW methods in conjunction with the Yee et al. and Yee & Sjögreen nonlinear filter approach.

The relative stability and accuracy are comparable among the KGP-type, DS-type split, ES, and ESSW methods. In certain test cases, ES and ESSW are more stable for a longer time integration. The KGP-type and DS-type methods are the least CPU-intensive for both the gas dynamics and MHD test cases. However, they do not conserve entropy. For the difficult Brio–Wu test case, unlike most published work reported in the literature, the current nonlinear filter approach, in conjunction of the various splitting methods, indicated that stability and high accuracy were obtained using very high-order spatial discretization and the Roe-type approximate Riemann solver. The numerical instability and highly spurious oscillatory solutions via standard high-order shock-capturing methods are drastically reduced using the ES and ESSW methods. In current studies, and according to our previous studies [19,20,22,109], EC, ECLOG, and ECLOGKP are the most CPU-intensive methods among the nine methods. They use at least twice the CPU per time step as the ES, ESDS, and ESSW methods; yet, they exhibit similar resolutions. DS is the least CPU-intensive.

8. Generalization to a Wider Class of Entropy Split Methods for Compressible Ideal MHD [23]

Since this is a new development, the majority of the content from [23] is repeated here. The new method consists of a two-point numerical flux portion and a non-conservative portion in such a way that entropy conservation holds without requiring the homogeneity property of the compressible inviscid flux. For high-order classical spatial central, DRP, or Padé (compact) spatial discretizations, this new approach can be proven to be entropy-conservative while, at the same time, allowing a wider class of symmetrizable inviscid flux derivatives. More importantly, we extend this new approach to derive an entropy split method for the equations of ideal MHD (which does not have homogeneous fluxes) using the Godunov symmetrizable ideal MHD formulation.

First, we present a new formulation of the entropy split method that is entropy-conserving but that relaxes the flux homogeneity requirement. Then, we derive an entropy split method that is entropy-conserving for the equations of MHD.

8.1. New Approach to High-Order Entropy Split Methods Using Two-Point Numerical Flux Approach

Consider the split form of a second-order, accurate, semi-discrete conservation law,

$$\Delta x\frac{d{\mathbf{q}}_j}{dt}+\omega ({\mathbf{h}}_{j+1/2}-{\mathbf{h}}_{j-1/2})+(1-\omega )A_j\frac{{\mathbf{v}}_{j+1}-{\mathbf{v}}_{j-1}}{2}=\mathbf{0},j=1,\dots ,N$$

(82)

where ${\mathbf{h}}_{j+1/2}$ is a numerical flux function consistent up to the second order with the flux $\mathbf{f}$ of the continuous problem, and $A=\partial \mathbf{f}/\partial \mathbf{v}$. The constant weight $\omega$, with $0\le \omega \le 1$, is the fraction of the conservative part of the discretization. We achieve the following result.

Theorem 4.

Assume that the numerical fluxes ${\mathbf{h}}_{j+1/2}$ of (82) are determined such that

$$\omega {({\mathbf{v}}_{j+1}-{\mathbf{v}}_j)}^T{\mathbf{h}}_{j+1/2}=(1-\omega )\frac{1}{2}({\mathbf{v}}_{j+1}^T{A}_{j+1}+{\mathbf{v}}_j^T{A}_j)({\mathbf{v}}_{j+1}-{\mathbf{v}}_j)+{\psi }_{j+1}-{\psi }_j,$$

(83)

where $\psi =\omega {\mathbf{v}}^T\mathbf{f}-F$. Then, the discretization (82) implies the entropy conservation law

$$\Delta x\frac{d}{dt}{E}_j\left(t\right)+H_{j+1/2}-H_{j-1/2}=0,$$

(84)

where the discrete entropy fluxes are given via

$$H_{j+1/2}=\frac{\omega }{2}{({\mathbf{v}}_{j+1}+{\mathbf{v}}_j)}^T{\mathbf{h}}_{j+1/2}-\frac{1-\omega }{4}({\mathbf{v}}_{j+1}^T{A}_{j+1}-{\mathbf{v}}_j^T{A}_j)({\mathbf{v}}_{j+1}-{\mathbf{v}}_j)-\frac{1}{2}({\psi }_{j+1}+{\psi }_j).$$

(85)

Proof. 

We multiply (82) by ${\mathbf{v}}_j^T$. This gives

$$\Delta x\frac{d{E}_j}{dt}+\omega {\mathbf{v}}_j^T({\mathbf{h}}_{j+1/2}-{\mathbf{h}}_{j-1/2})+(1-\omega ){\mathbf{v}}_j^T{A}_j\frac{{\mathbf{v}}_{j+1}-{\mathbf{v}}_{j-1}}{2}=\mathbf{0}.$$

The flux difference term is rewritten:

$$\begin{array}{c}{\mathbf{v}}_j^T({\mathbf{h}}_{j+1/2}-{\mathbf{h}}_{j-1/2})=\frac{1}{2}{({\mathbf{v}}_{j+1}+{\mathbf{v}}_j)}^T{\mathbf{h}}_{j+1/2}-\frac{1}{2}{({\mathbf{v}}_j+{\mathbf{v}}_{j-1})}^T{\mathbf{h}}_{j-1/2}-\hfill \\ \hfill \frac{1}{2}{({\mathbf{v}}_{j+1}-{\mathbf{v}}_j)}^T{\mathbf{h}}_{j+1/2}-\frac{1}{2}{({\mathbf{v}}_j-{\mathbf{v}}_{j-1})}^T{\mathbf{h}}_{j-1/2}.\end{array}$$

(86)

The non-conservative term is rewritten:

$$\begin{array}{c}{\mathbf{v}}_j^T{A}_j({\mathbf{v}}_{j+1}-{\mathbf{v}}_{j-1})=\hfill \\ \frac{1}{2}({\mathbf{v}}_{j+1}^T{A}_{j+1}+{\mathbf{v}}_j^T{A}_j)({\mathbf{v}}_{j+1}-{\mathbf{v}}_j)+\frac{1}{2}({\mathbf{v}}_j^T{A}_j+{\mathbf{v}}_{j-1}^T{A}_{j-1})({\mathbf{v}}_j-{\mathbf{v}}_{j-1})-\\ \frac{1}{2}({\mathbf{v}}_{j+1}^T{A}_{j+1}-{\mathbf{v}}_j^T{A}_j)({\mathbf{v}}_{j+1}-{\mathbf{v}}_j)+\frac{1}{2}({\mathbf{v}}_j^T{A}_j-{\mathbf{v}}_{j-1}^T{A}_{j-1})({\mathbf{v}}_j-{\mathbf{v}}_{j-1}).\end{array}$$

(87)

Adding together and collecting terms that are different gives

$$\begin{array}{c}\Delta x\frac{d{E}_j}{dt}+{\Delta }_{+}(\omega \frac{1}{2}{({\mathbf{v}}_j+{\mathbf{v}}_{j-1})}^T{\mathbf{h}}_{j-1/2}-(1-\omega )\frac{1}{4}({\mathbf{v}}_j^T{A}_j-{\mathbf{v}}_{j-1}^T{A}_{j-1})({\mathbf{v}}_j-{\mathbf{v}}_{j-1}))\hfill \\ -\omega \frac{1}{2}{({\mathbf{v}}_{j+1}-{\mathbf{v}}_j)}^T{\mathbf{h}}_{j+1/2}+(1-\omega )\frac{1}{4}({\mathbf{v}}_{j+1}^T{A}_{j+1}+{\mathbf{v}}_j^T{A}_j)({\mathbf{v}}_{j+1}-{\mathbf{v}}_j)\\ \hfill -\omega \frac{1}{2}{({\mathbf{v}}_j-{\mathbf{v}}_{j-1})}^T{\mathbf{h}}_{j-1/2}+(1-\omega )\frac{1}{4}({\mathbf{v}}_j^T{A}_j+{\mathbf{v}}_{j-1}^T{A}_{j-1})({\mathbf{v}}_j-{\mathbf{v}}_{j-1})=0,\end{array}$$

(88)

where the first line is of a conservative form. The forward difference is defined by ${\Delta }_{+}\left(u_j\right)=u_{j+1}-u_j$. The second and third lines can, by (83), be replaced with $-({\psi }_{j+1}-{\psi }_j)/2$ and $-({\psi }_j-{\psi }_{j-1})/2$, respectively. We obtain

$$\begin{array}{c}\Delta x\frac{d{E}_j}{dt}+{\Delta }_{+}(\omega \frac{1}{2}{({\mathbf{v}}_j+{\mathbf{v}}_{j-1})}^T{\mathbf{h}}_{j-1/2}-(1-\omega )\frac{1}{4}({\mathbf{v}}_j^T{A}_j-{\mathbf{v}}_{j-1}^T{A}_{j-1})({\mathbf{v}}_j-{\mathbf{v}}_{j-1}))\hfill \\ \hfill -\frac{1}{2}({\psi }_{j+1}-{\psi }_j)-\frac{1}{2}({\psi }_j-{\psi }_{j-1})=0.\end{array}$$

(89)

Entropy conservation with the numerical entropy flux (85) follows from (89) by rewriting the $\psi$-terms as

$$-\frac{1}{2}({\psi }_{j+1}+{\psi }_j)+\frac{1}{2}({\psi }_j+{\psi }_{j-1})=-\frac{1}{2}{\Delta }_{+}({\psi }_j+{\psi }_{j-1}).$$

This shows (84) and thereby proves the theorem. □

According to Theorem 4, all that is needed in order to define an entropy split method of form (82) that conserves entropy is to determine the numerical flux function, ${\mathbf{h}}_{j+1/2}$, so that it satisfies (83). As an example, this is next carried out for the equations of gas dynamics, using the flux (27) and the entropy (29). For the thermally perfect gas dynamics equations, the entropy split form (78) is known to conserve entropy. This example will demonstrate that (78) can be derived using Theorem 4.

The entropy variables corresponding to the Harten entropy functions (29) are

$$\mathbf{v}=\frac{\rho }{p}{s}^{\frac{1}{\alpha +\gamma }}(-\frac{\alpha }{\gamma -1}\frac{p}{\rho }-\frac{1}{2}{\left|\mathbf{u}\right|}^2{,u,v,w,-1)}^T.$$

(90)

Multiplying together (90) and the inviscid flux gives

$${\mathbf{v}}^T\mathbf{f}=-\frac{\alpha +1}{\gamma -1}\rho \widehat{u}{s}^{\frac{1}{\alpha +\gamma }},$$

so that

$$\psi =\omega {\mathbf{v}}^T\mathbf{f}-F=-\omega \frac{\alpha +1}{\gamma -1}\rho \widehat{u}{s}^{\frac{1}{\alpha +\gamma }}-\widehat{u}E=\omega \frac{\alpha +1}{\alpha +\gamma }\widehat{u}E-\widehat{u}E,$$

where the second equality follows from the entropy definition (29). If we take $\omega =\frac{\alpha +\gamma }{\alpha +1}=\beta /(\beta +1)$, then $\psi =0$. Furthermore, by using $\psi =0$ and ${\mathbf{v}}^{T}A=\beta \mathbf{f}$ (which can be shown by directly evaluating the product ${\mathbf{v}}^{T}A$) (83), we obtain the following:

$${({\mathbf{v}}_{j+1}-{\mathbf{v}}_j)}^T{\mathbf{h}}_{j+1/2}=\frac{1}{2}{({\mathbf{v}}_{j+1}-{\mathbf{v}}_j)}^T({\mathbf{f}}_{j+1}+{\mathbf{f}}_j),$$

which, clearly, can be satisfied by defining ${\mathbf{h}}_{j+1/2}$ as the standard second-order accurate centered flux, i.e., ${\mathbf{h}}_{j+1/2}=({\mathbf{f}}_j+{\mathbf{f}}_{j+1})/2$. In this sense, the entropy split method follows from (83) without explicit use of homogeneity. A by-product of this way of deriving the entropy split method is that, thanks to Theorem 4, the method is consistent with the local entropy conservation law (84) with a numerical entropy flux

$$H_{j+1/2}=\frac{\beta }{(\beta +1)}\frac{1}{4}({({\mathbf{v}}_{j+1}+{\mathbf{v}}_j)}^T({\mathbf{f}}_{j+1}+{\mathbf{f}}_j)-{({\mathbf{v}}_{j+1}-{\mathbf{v}}_j)}^T({\mathbf{f}}_{j+1}-{\mathbf{f}}_j)).$$

This local entropy conservation is not obtained from the standard SBP analysis leading to (34).

Hence, for fluxes that are not homogeneous, it would be possible to derive an entropy split method by determining numerical fluxes that satisfy (83). Such a method would be more complicated since the conservative part of the split formula would not necessarily be a standard centered approximation. This is a subject for a future investigation.

The above shows that the proofs used the second-order classical spatial center methods. In the same procedure as indicated in Sjögreen & Yee and Sjögreen et al. [18,19,20,21], it is straightforward to extend the new formulation to higher than second-order classical central, DRP, and Padé spatial discretizations. We next investigate the extension the two-point numerical flux entropy split method that is entropy-conserving to the equations of ideal MHD with fluxes that do not have the homogeneity property.

8.2. Extension of the New Entropy Split Formulation for the Equations of Ideal MHD

The fluxes of the equations of ideal MHD are not homogeneous. Moreover, the flux Jacobians with respect to the entropy variables are not symmetric. This means that the standard technique used for gas dynamics to prove that (78) conserves entropy does not directly carry over to the corresponding entropy split form of the equations of MHD. This section uses an extension of Theorem 4 to the equations of ideal MHD to define an entropy-conserving entropy split approximation. Unlike the gas dynamics case, the numerical flux of the conservative part, ${\mathbf{h}}_{j+1/2}$, of the MHD entropy split approximation is not the standard centered flux.

Consider the ideal MHD system,

$${\mathbf{q}}_t+{\mathbf{f}}_x+{\mathbf{g}}_y+{\mathbf{h}}_z+\mathbf{e}\left(\nabla ·\mathbf{B}\right)=\mathbf{0}$$

with conserved variables

$$\mathbf{q}={(\rho ,\rho u_1,\rho u_2,\rho u_3,e,B_1,B_2,B_3)}^T,$$

and introduce the notation $\mathbf{u}={\left(u_1{u}_2{u}_3\right)}^T$, $\mathbf{B}={\left(B_1{B}_2{B}_3\right)}^T$, ${\left|\mathbf{u}\right|}^2=u_1^2+u_2^2+u_3^2$, and ${\left|\mathbf{B}\right|}^2=B_1^2+B_2^2+B_3^2$. The flux in direction $(k_1,k_2,k_3)$ is

$$\mathbf{f}=\left(\begin{array}{c}\rho \widehat{u}\\ \rho \widehat{u}{u}_1+k_1(p+\frac{1}{2}{\left|\mathbf{B}\right|}^2)-\widehat{B}{B}_1\\ \rho \widehat{u}{u}_2+k_2(p+\frac{1}{2}{\left|\mathbf{B}\right|}^2)-\widehat{B}{B}_2\\ \rho \widehat{u}{u}_3+k_3(p+\frac{1}{2}{\left|\mathbf{B}\right|}^2)-\widehat{B}{B}_3\\ \widehat{u}(e+p+\frac{1}{2}{\left|\mathbf{B}\right|}^2)-\widehat{B}{\mathbf{u}}^T\mathbf{B}\\ \widehat{u}{B}_1-\widehat{B}{u}_1\\ \widehat{u}{B}_2-\widehat{B}{u}_2\\ \widehat{u}{B}_3-\widehat{B}{u}_3\end{array}\right),$$

(91)

where $\widehat{u}=k_1{u}_1+k_2{u}_2+k_3{u}_3$ and $\widehat{B}=k_1{B}_1+k_2{B}_2+k_3{B}_3$. The “source term” vector is

$$\mathbf{e}={(0,B_1,B_2,B_3,{\mathbf{u}}^T\mathbf{B},u_1,u_2,u_3)}^T,$$

which was used in [118]. Other choices for $\mathbf{e}$ are possible. The fluxes do not have the homogeneity property. The total energy is related to the pressure, p, via the constitutive law,

$$e=\frac{p}{\gamma -1}+\frac{1}{2}{\rho \left|\mathbf{u}\right|}^2+\frac{1}{2}{\left|\mathbf{B}\right|}^2.$$

Here, we use the Harten entropies (29). It is straightforward to verify that $E_{\mathbf{q},\mathbf{q}}>0$ for $\alpha <-\gamma$ or $\alpha >0$ for the equations of MHD. Similarly, the entropy flux is $F=\widehat{u}E$, just as in the gas dynamics case. The entropy variables $\mathbf{v}=E_{\mathbf{q}}$ are found via differentiation to be

$$\mathbf{v}=\frac{\rho }{p}{s}^{\frac{1}{\alpha +\gamma }}(-\frac{\alpha }{\gamma -1}\frac{p}{\rho }-\frac{1}{2}{\left|\mathbf{u}\right|}^2,u_1,u_2,u_3,-1,B_1,B_2,B_3{)}^T.$$

(92)

Consider the one-dimensional split discretization of the MHD equations,

$$\Delta x\frac{d{\mathbf{q}}_j}{dt}+\omega ({\mathbf{h}}_{j+1/2}-{\mathbf{h}}_{j-1/2})+(1-\omega )A_j\frac{{\mathbf{v}}_{j+1}-{\mathbf{v}}_{j-1}}{2}+\omega {\mathbf{e}}_j\frac{{\widehat{B}}_{j+1}-{\widehat{B}}_{j-1}}{2}=\mathbf{0},$$

(93)

for $j=1,\dots ,N$, where ${\mathbf{h}}_{j+1/2}$ is a numerical flux to be determined to make (93) conserve entropy. The symmetric matrix A is given via

$$A=\frac{\partial \mathbf{f}}{\partial \mathbf{v}}+C,$$

where C is the matrix obtained by rewriting $\mathbf{e}{\widehat{B}}_x=C{\mathbf{v}}_x$. A simple modification of Theorem 4 to include the MHD “source term” gives

Theorem 5.

Assume that the numerical fluxes ${\mathbf{h}}_{j+1/2}$ of (93) satisfy

$$\begin{array}{c}\omega {({\mathbf{v}}_{j+1}-{\mathbf{v}}_j)}^T{\mathbf{h}}_{j+1/2}=\frac{1-\omega }{2}({\mathbf{v}}_{j+1}^T{A}_{j+1}+{\mathbf{v}}_j^T{A}_j)({\mathbf{v}}_{j+1}-{\mathbf{v}}_j)\hfill \\ +\frac{\omega }{2}({\mathbf{v}}_{j+1}^T{\mathbf{e}}_{j+1}+{\mathbf{v}}_j^T{\mathbf{e}}_j)({\widehat{B}}_{j+1}-{\widehat{B}}_j)+{\psi }_{j+1}-{\psi }_j,\end{array}$$

(94)

where $\psi =\omega {\mathbf{v}}^T\mathbf{f}-F$. Then, the discretization (93) implies the entropy conservation law

$$\Delta x\frac{E_j}{dt}+H_{j+1/2}-H_{j-1/2}=0,$$

where the discrete entropy fluxes are given via

$$\begin{array}{c}{H}_{j+1/2}=\frac{\omega }{2}{({\mathbf{v}}_{j+1}+{\mathbf{v}}_j)}^T{\mathbf{h}}_{j+1/2}-\frac{1-\omega }{4}({\mathbf{v}}_{j+1}^T{A}_{j+1}-{\mathbf{v}}_j^T{A}_j)({\mathbf{v}}_{j+1}-{\mathbf{v}}_j)\hfill \\ \hfill -\frac{\omega }{4}({\mathbf{v}}_{j+1}^T{\mathbf{e}}_{j+1}-{\mathbf{v}}_j^T{\mathbf{e}}_j)({\widehat{B}}_{j+1}-{\widehat{B}}_j)-\frac{1}{2}({\psi }_{j+1}+{\psi }_j).\end{array}$$

(95)

To see what this means for the equations of MHD, it is straightforward to evaluate

$${\mathbf{v}}^T\mathbf{f}=-\frac{\alpha +1}{\gamma -1}\rho \widehat{u}{s}^{\frac{1}{\alpha +\gamma }}+\frac{\rho }{p}{s}^{\frac{1}{\alpha +\gamma }}\left(\frac{1}{2}\widehat{u}{\left|\mathbf{B}\right|}^2-{\mathbf{u}}^T\mathbf{B}\widehat{B}\right).$$

The weighted entropy flux potential becomes

$$\psi =\omega {\mathbf{v}}^T\mathbf{f}-F=(-\omega \frac{\alpha +1}{\gamma -1}+\frac{\alpha +\gamma }{\gamma -1})\rho \widehat{u}{s}^{\frac{1}{\alpha +\gamma }}+\omega \frac{\rho }{p}{s}^{\frac{1}{\alpha +\gamma }}\left(\frac{1}{2}\widehat{u}{\left|\mathbf{B}\right|}^2-\widehat{B}{\mathbf{u}}^T\mathbf{B}\right).$$

Setting $\omega =(\alpha +\gamma )/(\alpha +1)=\beta /(\beta +1)$ eliminates the first term of the entropy flux potential. We denote

$${\psi }^M=\frac{\rho }{p}{s}^{\frac{1}{\alpha +\gamma }}\left(\frac{1}{2}\widehat{u}{\left|\mathbf{B}\right|}^2-\widehat{B}{\mathbf{u}}^T\mathbf{B}\right).$$

(96)

We split the MHD flux (91) into a gas dynamics part and a magnetic part, $\mathbf{f}={\mathbf{f}}^G+{\mathbf{f}}^M$, with

$${\mathbf{f}}^G=\left(\begin{array}{c}\rho \widehat{u}\\ \rho \widehat{u}{u}_1+k_{1}p\\ \rho \widehat{u}{u}_2+k_{2}p\\ \rho \widehat{u}{u}_3+k_{3}p\\ \widehat{u}\left(e+p-\frac{1}{2}{\left|\mathbf{B}\right|}^2\right)\\ 0\\ 0\\ 0\end{array}\right){\mathbf{f}}^M=\left(\begin{array}{c}0\\ k_1\frac{1}{2}{\left|\mathbf{B}\right|}^2-\widehat{B}{B}_1\\ k_2\frac{1}{2}{\left|\mathbf{B}\right|}^2-\widehat{B}{B}_2\\ k_3\frac{1}{2}{\left|\mathbf{B}\right|}^2-\widehat{B}{B}_3\\ \widehat{u}{\left|\mathbf{B}\right|}^2-\widehat{B}{\mathbf{u}}^T\mathbf{B}\\ \widehat{u}{B}_1-\widehat{B}{u}_1\\ \widehat{u}{B}_2-\widehat{B}{u}_2\\ \widehat{u}{B}_3-\widehat{B}{u}_3\end{array}\right).$$

Then, the following holds.

Theorem 6.

Let the numerical flux of the entropy split method (93) be selected as

$${\mathbf{h}}_{j+1/2}=\frac{1}{2}({\mathbf{f}}_{j+1}^G+{\mathbf{f}}_j^G)+{\mathbf{h}}_{j+1/2}^M,$$

(97)

where the magnetic part of the numerical flux, ${\mathbf{h}}_{j+1/2}^M$, satisfies

$${({\mathbf{v}}_{j+1}-{\mathbf{v}}_j)}^T{\mathbf{h}}_{j+1/2}^M=\frac{1}{2}({\mathbf{v}}_{j+1}^T{\mathbf{e}}_{j+1}+{\mathbf{v}}_j^T{\mathbf{e}}_j)({\widehat{B}}_{j+1}-{\widehat{B}}_j)+{\psi }_{j+1}^M-{\psi }_j^M.$$

(98)

Then, (94) is satisfied. Hence, according to Theorem 5 the discretization (93) with the numerical flux (97) is an entropy-conserving discretization.

Proof. 

To apply (94), the quantity ${\mathbf{v}}^{T}A$ is required. A straightforward but lengthy evaluation gives

$$\begin{array}{c}{\mathbf{v}}^{T}A=-\frac{\alpha +\gamma }{\gamma -1}(\rho \widehat{u},\rho \widehat{u}{u}_1+k_{1}p,\rho \widehat{u}{u}_2+k_{2}p,\rho \widehat{u}{u}_3+k_{3}p,\widehat{u}(e+p-\frac{1}{2}{\left|\mathbf{B}\right|}^2),0,0,0)\hfill \\ \hfill =\beta {\mathbf{f}}^G,\end{array}$$

(99)

because $(\alpha +\gamma )/(\alpha +1)=\beta /(\beta +1)$ implies $(\alpha +\gamma )/(\gamma -1)=-\beta$. Dividing (94) by $\omega$, using $(1-\omega )\beta =\omega$ and (99), allows us to rewrite (94) as

$$\begin{array}{c}{({\mathbf{v}}_{j+1}-{\mathbf{v}}_j)}^T{\mathbf{h}}_{j+1/2}=\frac{1}{2}({\mathbf{f}}_{j+1}^G+{\mathbf{f}}_j^G)({\mathbf{v}}_{j+1}-{\mathbf{v}}_j)+\frac{1}{2}({\mathbf{v}}_{j+1}^T{\mathbf{e}}_{j+1}+{\mathbf{v}}_j^T{\mathbf{e}}_j)({\widehat{B}}_{j+1}-{\widehat{B}}_j)\hfill \\ \hfill +{\psi }_{j+1}^M-{\psi }_j^M.\end{array}$$

(100)

Inserting the numerical flux (97) into the left-hand side of (100) gives (98), which shows that (94) is equivalent to (98) when (97) is used. □

8.3. Derivation of a Split-Form Entropy-Conserving Method

An explicit expression for ${\mathbf{h}}_{j+1/2}^M$ that satisfies (98) is derived by introducing the parameter vector

$$\mathbf{z}=(\frac{\rho }{p}{s}^{\frac{1}{\alpha +\gamma }},u_1,u_2,u_3,p,B_1,B_2,B_3).$$

The entropy variables and entropy flux potential are written in terms of the parameter vector as

$$\begin{array}{cc}\hfill v_1& =-\frac{1}{2}{z}_1(z_2^2+z_3^2+z_4^2)-\frac{\alpha }{\gamma -1}{z}_1^{-\frac{\gamma }{\alpha }}{z}_5^{-\frac{\gamma -1}{\alpha }}\hfill \\ \hfill v_2& =z_1{z}_2{v}_3=z_1{z}_3{v}_4=z_1{z}_4\hfill \\ \hfill v_5& =-z_1\hfill \\ \hfill v_6& =z_1{z}_6{v}_7=z_1{z}_7{v}_8=z_1{z}_8\hfill \\ \hfill {\psi }^M& =\frac{1}{2}\widehat{z}{z}_1(z_6^2+z_7^2+z_8^2)-\tilde{z}{z}_1(z_2{z}_6+z_3{z}_7+z_4{z}_8)\hfill \end{array}$$

where $\widehat{z}=k_1{z}_2+k_2{z}_3+k_3{z}_4$ and $\tilde{z}=k_1{z}_6+k_2{z}_7+k_3{z}_8$. Furthermore,

$${\mathbf{e}}^T\mathbf{v}=z_1(z_2{z}_6+z_3{z}_7+z_4{z}_8).$$

We determine a magnetic numerical flux function, ${\mathbf{h}}_{j+1/2}^M=\mathbf{h}({\mathbf{u}}_{j+1},{\mathbf{u}}_{\mathbf{j}})$, so that it satisfies (98), expressed as

$${\left({\mathbf{h}}_{j+1/2}^M\right)}^T\Delta \mathbf{v}=\left\{{\mathbf{v}}^T\mathbf{e}\right\}\Delta \widehat{B}+\Delta {\psi }^M,$$

(101)

where the difference and average are denoted as

$$\Delta u=(u_{j+1}-u_j)\mathrm{and}\left\{u\right\}=(u_{j+1}+u_j)/2,$$

respectively.

Expressed in the intermediate variables, the right-hand side of (101) becomes

$$\left\{z_1(z_2{z}_6+z_3{z}_7+z_4{z}_8)\right\}\Delta \tilde{z}+\Delta \left(\frac{1}{2}\widehat{z}{z}_1(z_6^2+z_7^2+z_8^2)\right)-\Delta \left(\tilde{z}{z}_1(z_2{z}_6+z_3{z}_7+z_4{z}_8)\right),$$

which we expand as

$$\begin{array}{c}\frac{1}{2}\left\{\widehat{z}\right\}\{z_6^2+z_7^2+z_8^2\}\Delta z_1+\left\{\widehat{z}{z}_1\right\}(\left\{z_6\right\}\Delta z_6+\left\{z_7\right\}\Delta z_7+\left\{z_8\right\}\Delta z_8)\hfill \\ +\frac{1}{2}\left\{z_1\right\}\{z_6^2+z_7^2+z_8^2\}\Delta \widehat{z}-\left\{\tilde{z}\right\}(\left\{z_2\right\}\left\{z_6\right\}+\left\{z_3\right\}\left\{z_7\right\}+\left\{z_4\right\}\left\{z_8\right\})\Delta z_1\\ -\left\{\tilde{z}\right\}\left\{z_1\right\}(\left\{z_6\right\}\Delta z_2+\left\{z_7\right\}\Delta z_3+\left\{z_8\right\}\Delta z_4)\\ \hfill -\left\{\tilde{z}\right\}(\left\{z_1{z}_2\right\}\Delta z_6+\left\{z_1{z}_3\right\}\Delta z_7+\left\{z_1{z}_4\right\}\Delta z_8).\end{array}$$

(102)

For the left-hand side of (101), the entropy variables are written out in terms of $\mathbf{z}$ to obtain

$$\begin{array}{c}{h}^{\left(2\right)}\left\{z_1\right\}\Delta z_2+h^{\left(3\right)}\left\{z_1\right\}\Delta z_3+h^{\left(4\right)}\left\{z_1\right\}\Delta z_4\hfill \\ +\left(-h^{\left(5\right)}+\left\{z_2\right\}{h}^{\left(2\right)}+\left\{z_3\right\}{h}^{\left(3\right)}+\left\{z_4\right\}{h}^{\left(4\right)}+\left\{z_6\right\}{h}^{\left(6\right)}+\left\{z_7\right\}{h}^{\left(7\right)}+\left\{z_8\right\}{h}^{\left(8\right)}\right)\Delta z_1\\ \hfill +h^{\left(6\right)}\left\{z_1\right\}\Delta z_6+h^{\left(7\right)}\left\{z_1\right\}\Delta z_7+h^{\left(8\right)}\left\{z_1\right\}\Delta z_8,\end{array}$$

(103)

where $h^{\left(j\right)}$, $j=1,\dots ,8$ denotes the components of the magnetic numerical flux, ${\mathbf{h}}^M$. By equating $\Delta z_j$ terms, the magnetic part of the numerical flux is found to be

$${\mathbf{h}}^M=\left(\begin{array}{c}0\\ \frac{1}{2}{k}_1\{z_6^2+z_7^2+z_8^2\}-\left\{\tilde{z}\right\}\left\{z_6\right\}\\ \frac{1}{2}{k}_2\{z_6^2+z_7^2+z_8^2\}-\left\{\tilde{z}\right\}\left\{z_7\right\}\\ \frac{1}{2}{k}_3\{z_6^2+z_7^2+z_8^2\}-\left\{\tilde{z}\right\}\left\{z_8\right\}\\ \frac{\left\{z_1\widehat{z}\right\}}{\left\{z_1\right\}}({\left\{z_6\right\}}^2+{\left\{z_7\right\}}^2+{\left\{z_8\right\}}^2)-\left\{\tilde{z}\right\}(\frac{\left\{z_1{z}_2\right\}}{\left\{z_1\right\}}\left\{z_6\right\}+\frac{\left\{z_1{z}_3\right\}}{\left\{z_1\right\}}\left\{z_7\right\}+\frac{\left\{z_1{z}_4\right\}}{\left\{z_1\right\}}\left\{z_8\right\})\\ \frac{\left\{z_1\widehat{z}\right\}}{\left\{z_1\right\}}\left\{z_6\right\}-\left\{\tilde{z}\right\}\frac{\left\{z_1{z}_2\right\}}{\left\{z_1\right\}}\\ \frac{\left\{z_1\widehat{z}\right\}}{\left\{z_1\right\}}\left\{z_7\right\}-\left\{\tilde{z}\right\}\frac{\left\{z_1{z}_3\right\}}{\left\{z_1\right\}}\\ \frac{\left\{z_1\widehat{z}\right\}}{\left\{z_1\right\}}\left\{z_8\right\}-\left\{\tilde{z}\right\}\frac{\left\{z_1{z}_4\right\}}{\left\{z_1\right\}}\end{array}\right).$$

In terms of the original variables, the flux is

$${\mathbf{h}}^M=\left(\begin{array}{c}0\\ \frac{1}{2}{k}_1\{B_1^2+B_2^2+B_3^2\}-\left\{\widehat{B}\right\}\left\{B_1\right\}\\ \frac{1}{2}{k}_2\{B_1^2+B_2^2+B_3^2\}-\left\{\widehat{B}\right\}\left\{B_2\right\}\\ \frac{1}{2}{k}_3\{B_1^2+B_2^2+B_3^2\}-\left\{\widehat{B}\right\}\left\{B_3\right\}\\ \frac{\left\{z_1\widehat{u}\right\}}{\left\{z_1\right\}}({\left\{B_1\right\}}^2+{\left\{B_2\right\}}^2+{\left\{B_3\right\}}^2)-\left\{\widehat{B}\right\}(\frac{\left\{z_1{u}_1\right\}}{\left\{z_1\right\}}\left\{B_1\right\}+\frac{\left\{z_1{u}_2\right\}}{\left\{z_1\right\}}\left\{B_2\right\}+\frac{\left\{z_1{u}_3\right\}}{\left\{z_1\right\}}\left\{B_3\right\})\\ \frac{\left\{z_1\widehat{u}\right\}}{\left\{z_1\right\}}\left\{B_1\right\}-\left\{\widehat{B}\right\}\frac{\left\{z_1{u}_1\right\}}{\left\{z_1\right\}}\\ \frac{\left\{z_1\widehat{u}\right\}}{\left\{z_1\right\}}\left\{B_2\right\}-\left\{\widehat{B}\right\}\frac{\left\{z_1{u}_2\right\}}{\left\{z_1\right\}}\\ \frac{\left\{z_1\widehat{u}\right\}}{\left\{z_1\right\}}\left\{B_3\right\}-\left\{\widehat{B}\right\}\frac{\left\{z_1{u}_3\right\}}{\left\{z_1\right\}}\end{array}\right),$$

with $z_1=\frac{\rho }{p}{s}^{\frac{1}{\alpha +\gamma }}$.

Hence, the entropy split scheme for MHD (93) with a numerical flux function,

$$\mathbf{h}=\left\{{\mathbf{f}}^G\right\}+{\mathbf{h}}^M,$$

(104)

and with $\omega =\beta /(\beta +1)$ has mathematically strict entropy conservation. Note that only the gas dynamics part, ${\mathbf{f}}^G$, of the numerical flux (104) is a standard-centered flux. The magnetic field part of (104) is of a form similar to a Tadmor-type entropy-conserving numerical flux.

Using the two-point flux (104) in (93) leads to a second-order accurate method. A generalization to a higher order is done in the standard way [110]. The pth-order accurate generalization of (93), for $p=2s$ even, explicitly written out, is

$$\begin{array}{c}\frac{d{\mathbf{q}}_j}{dt}+\frac{\beta }{\beta +1}\left(D^{\left(p\right)}{\mathbf{f}}_j^G+\frac{1}{\Delta x}\sum _{k=1}^{s}2{\alpha }_k^{\left(p\right)}({\mathbf{h}}^M({\mathbf{q}}_{j+k},{\mathbf{q}}_j)-{\mathbf{h}}^M({\mathbf{q}}_j,{\mathbf{q}}_{j-k}))\right)\hfill \\ +\frac{1}{\beta +1}{A}_j{D}^{\left(p\right)}{\mathbf{v}}_j+\frac{\beta }{\beta +1}{\mathbf{e}}_j{D}^{\left(p\right)}{\widehat{B}}_j=\mathbf{0},\end{array}$$

(105)

where

$$D^{\left(p\right)}{u}_j=\frac{1}{\Delta x}\sum _{k=1}^s{\alpha }_k^{\left(p\right)}(u_{j+k}-u_{j-k})$$

denotes the pth-order accurate centered difference operator. The elements of the symmetric Jacobian matrix $A_j$ are given in Appendix A of [30].

When the magnetic field vanishes, the equations of ideal MHD become the inviscid Euler equations of gas dynamics. For thermally perfect gas, the new entropy split scheme in this case becomes the standard entropy split scheme [26] obtained by setting the magnetic flux ${\mathbf{h}}_{j+1/2}^M=\mathbf{0}$ in (105). Entropy splitting for thermally perfect gas dynamics equations and extensions to deal with non-conservation were studied in [22]. For test cases with extensive grid refinement and eight different high-order methods’ comparison, including the long-time integration of shock-free turbulence and turbulence with shocks, see [22].

8.4. Numerical Experiments

This section demonstrates the performance of the new entropy split method, (93) with $\omega =\beta /(\beta +1)$ and numerical flux (104), in a few simple MHD examples. The new method is denoted as ESnew. Method (93) with the centered flux ${\mathbf{h}}_{j+1/2}=({\mathbf{f}}_{j+1}+{\mathbf{f}}_j)/2$ is used for comparison. This method, which was described and evaluated using numerical examples in [30], is denoted as ES. Note that ES does not conserve entropy, whereas ESnew does.

The spatial discretization used in the numerical examples is of order accuracy eight, given (in one space dimension) for ESnew via (105) with $p=8$. Unless otherwise stated, the semi-discrete equations are integrated in time with a fourth-order accurate Runge–Kutta method using cfl-number 0.4.

8.4.1. Alfvén Wave

This example tests the stability of the methods by translating an Alfvén wave for a long time. Eventually, nonlinear instabilities develop, and the solution accuracy breaks down. The Alfvén wave has a constant density and pressure, $\rho =p=1$. The velocity and magnetic field are given via

$$\begin{array}{cc}\hfill \mathbf{u}& =A(-sin\alpha sin\varphi ,cos\alpha sin\varphi ,cos\varphi )\hfill \end{array}$$

(106)

$$\begin{array}{cc}\hfill \mathbf{B}& =(cos\alpha -Asin\alpha sin\varphi ,sin\alpha +Acos\alpha sin\varphi ,Acos\varphi ),\hfill \end{array}$$

(107)

where the constant angle $\alpha ={30}^{\circ }$, the constant amplitude $A=0.1$, and the phase is $\varphi =2\pi (xcos\alpha +ysin\alpha +t)$. The heat capacity ratio is $\gamma =5/3$. The wave is computed on a two-dimensional domain of size 1.3546 × 2.1915, using periodic boundary conditions and grid spacing of $\Delta x=\Delta y=0.022$, corresponding to $62\times 104$ grid points. The equations are integrated to time 200, which corresponds to approximately 69,000 time steps.

Figure 32 shows the maximum norm of the error vs. time for the ESnew and ES schemes for different values of $\beta$. It is seen that ESnew keeps the error smaller for $\beta =2.5$ and $\beta =25.5$, but ES has a lower error when $\beta =1$. In the limit $\beta \to 0$, both methods are identical, being 100% non-conservative. As $\beta$ grows, the conservative fraction of the discretization increases. In the limit $\beta \to \infty$, the discretizations become 100% conservative. In the conservative limit, ESnew uses the numerical flux (93), whereas ES uses the standard centered flux. This might explain why the largest improvement in stability with ESnew occurs for the largest value of $\beta$.

8.4.2. Brio–Wu MHD Riemann Problem

We selected the same Brio–Wu test case as in Section 7 with the same initial data, as it is very challenging for standard shock-capturing methods. The initial state vectors are repeated here.

$$\begin{array}{cc}\hfill ({\rho }_L,u_L,p_L,B^{\left(x\right)},B^{\left(y\right)},B^{\left(z\right)})& =(1,0,1,0.75,1,0)\hfill \end{array}$$

(108)

$$\begin{array}{cc}\hfill ({\rho }_R,u_R,p_R,B^{\left(x\right)},B^{\left(y\right)},B^{\left(z\right)})& =(0.125,0,0.1,0.75,-1,0).\hfill \end{array}$$

(109)

Due to the complex shock wave interaction, the nonlinear filter version of the new formulation of the entropy split method is employed [26,27,67]. First, a post-processing shock-capturing filter is applied after each time step in order to suppress spurious oscillations via the non-dissipative entropy split scheme. The filter uses a wavelet-based adaptive flow sensor, together with the dissipation of the seventh-order WENO scheme [61,119]. See [17,27,67] for details.

Then, a shock detector is used to adaptively switch to a completely conservative approximation in the neighborhood of shock waves in order to avoid errors in shock speed due to non-conservative effects. To this end, a sensor variable, $s_j$, satisfying $0\le s_j\le 1$, is introduced in the discretization (93), giving

$$\begin{array}{cc}\Delta x\frac{d{\mathbf{q}}_j}{dt}+\frac{\beta +s_j}{\beta +1}\left({\mathbf{h}}_{j+1/2}-{\mathbf{h}}_{j-1/2}+{\mathbf{e}}_j\frac{{\widehat{B}}_{j+1}-{\widehat{B}}_{j-1}}{2}\right)+\hfill & \\ & \hfill \frac{1-s_j}{\beta +1}{A}_j\frac{{\mathbf{v}}_{j+1}-{\mathbf{v}}_{j-1}}{2}=\mathbf{0},\end{array}$$

(110)

so that the entropy split method is obtained whenever $s_j=0$, and the fully conservative flux difference is obtained for $s_j=1$. Hence, $s_j$ should be designed to be one at shock waves and zero otherwise. The computations presented here used

$$s_j=\frac{{(1-r_j)}^6}{{(1+r_j)}^6},r_j=\frac{{\left({\Delta }_{+}{p}_j\right)}^2+ϵ}{{\left({\Delta }_{-}{p}_j\right)}^2+ϵ},$$

where $ϵ$ is a small number preventing a division by zero. See [22] for details.

The problem is solved on the domain $0\le x\le 1$ with the initial state

$$(\rho ,p,B_1,B_2)=\left\{\begin{array}{cc}(1,1,0.75,1)\hfill & x<0.5\hfill \\ (0.125,0.1,0.75,-1)\hfill & x\ge 0.5\hfill \end{array}\right..$$

The velocities and $B_3$ are zero on both sides of $x=0.5$. The heat capacity ratio is $\gamma =5/3$. The solution consists of, from left to right, a fast rarefaction wave, a compound slow wave, contact discontinuity, a slow shock, and a fast rarefaction wave. The waves are outlined in Figure 33.

The computations used 801 grid points and were run to time $0.1$. Figure 34 shows the density computed via ESSWnew and ESSW for $\beta =1.5,2.5$, and $25.5$. Figure 35 and Figure 36 similarly show the $B_2$ magnetic field component and the pressure, respectively. The exact solution, plotted in black, was obtained by solving the Rankine–Hugoniot condition across discontinuities and using a highly accurate ODE solver across the rarefaction waves.

It can be seen that both ESSW and ESSWnew capture the wave structure well. The close-up plots of the solution reveal that the ESSWnew method appears to be somewhat more oscillatory, but overall, the differences are minor. The solutions for $\beta =25.5$, which puts the most weight on the conservative part, appear to conform best with the exact solution for both ESSWnew and ESSW.

8.5. Shock–Oscillations Interaction MHD Problem

The purpose of this problem is to study the conservative properties of the methods. This is of interest because the entropy split scheme is not in a conservation form. For comparison, the results of the Tadmor-type entropy conservative scheme developed in [19], here denoted as ECH, are also shown. ECH preserves entropy (29) and is in a conservation form.

This test problem is a magnetohydrodynamic version of the Shu–Osher one-dimensional shock–oscillations interaction problem for gas dynamics. The magnetic field components were added to the original problem to meet two criteria. First, the shock speed is approximately the same as in the gas dynamics version of the problem. Secondly, no part of the flow features generated via the interaction exits at the boundaries, so the boundary values are constant.

The computational domain is $-5<x<5$, and the initial data are given via

$$\left(\begin{array}{c}\rho \\ u_1\\ u_2\\ p\\ B_1\\ B_2\end{array}\right)=\left\{\begin{array}{cc}\left(\begin{array}{c}2.824717\\ 2.293010\\ -0.132520\\ 8.934660\\ 1\\ 0.670399\end{array}\right)\hfill & x\le -4\hfill \\ \left(\begin{array}{c}1+0.2sin5x\\ 0\\ 0\\ 1\\ 1\\ 0.2\end{array}\right)\hfill & x>-4\hfill \end{array}\right.,$$

with $w=B_3=0$ on both sides. The problem is solved up to time 1.8.

Figure 37 shows the density, pressure, and y-component of the magnetic field at time 1.8 computed on a grid of 401 grid points using ESnew (red), ESSWnew (green), and ECH (blue). The reference solution (black) was obtained via the seventh-order accurate WENO scheme using 10,001 grid points. Figure 38 shows close-ups near the oscillatory region of the plots of Figure 37. Interestingly, the conservation errors of ESnew appear to lead to excess mass behind the shock wave, rather than an error in shock location. The conservatively modified ESSWnew appears to work well in correcting the conservation error.

Figure 39 compares the time evolution of the total mass for ESnew, ESSWnew, and ECH. The plotted total mass is compensated for the inflow and outflow across the domain boundaries and normalized to be zero initially, so the plotted quantity is $M^n-M^0$, where

$$M^n=\sum _{j=1}^N{\rho }_j^n\Delta x+t_n(f_R-f_L)$$

(111)

and where $f_R$ and $f_L$ are the constant mass fluxes across the right and left domain boundaries. It is straightforward to verify that $M^n$ is constant for numerical schemes of a conservative form. For plots of other conserved quantities, the density and mass fluxes in (111) are replaced with the appropriate conserved variable and flux.

Figure 39 clearly shows the non-conservativeness of ESnew. The introduction of ESSWnew significantly reduces the conservation error, and ECH shows perfect mass conservation, in agreement with the theory.

Execution Time Comparison

The computational efficiency of the schemes is strongly dependent on how they are implemented and the hardware used. We here give execution times using the solver ADPDIS3D, which provides a fortran implementation of the eight-order accurate version of the schemes. The code was compiled with the GNU fortran compiler and run on a single Intel i7-4770 CPU at 3.5 GHz. The CPU times reported were measured when running 3500 time steps of the Alfvén problem on the same grid as used in Section 8.4.1. The execution times in seconds for ESnew, ESSWnew, and ECH are presented in

Table 1. Execution times when running the Alfvén problem for 3500 time steps.

Method	CPU Seconds	Ratio
ECH	252	2.02
ESSWnew	127	1.02
ESnew	125	1

.

Table 1 shows that adding the shock sensor has a negligible effect on the execution time. Secondly, the Tadmor-type entropy conservative scheme has twice the computational cost of the entropy split method. A possible explanation for this is that, in the current implementation, ECH evaluates four exponentiations at each grid point, whereas the entropy split methods only use two exponentiations. The exponentiations are required when evaluating the entropy (29) and for ECH also to form an exponential average; see [30] for exact details of ECH.

9. Concluding Remark

This manuscript has presented a short historical overview and recent advances in high-order FDM development for DNS, LES, and ILES computations of compressible turbulence. The overview focused on recent trends in turbulence/shock numerical simulations that used the adaptive blending of high-order structure-preserving non-dissipative methods (classical central, Padé, and DRP) with high-order shock-capturing methods in such a way that high-order shock-capturing schemes are active only in the vicinity of shock/shear waves, high-gradient, and spurious high-frequency oscillation regions guided via flow sensors. Any efficient and high-resolution high-order shock-capturing methods can be part of the blending of method procedures. Typically, the adaptive blending of more than one method falls under two camps: hybrid methods and nonlinear filter methods. By examining the construction of these two approaches using the same non-dissipative high-order spatial discretization, the same high-order shock-capturing method, and the same flow sensor, the nonlinear filter approach was deemed to be more efficient and less CPU-intensive. These basic concepts of the two different approaches in the blending of more than one method are applicable to 3D FVM, FEM, DG, and spectral element methods.

The formulation of structure-preserving methods in a curvilinear grid was presented with a focus on the Yee et al. and Sjogreen & Yee Entropy Split Methods, which are entropy-conserving [19,23,26,30], since these entropy-conserving methods are not very familiar in the numerical method development community or among practitioners in the compressible turbulence simulation community.

According to the wide variety of test cases for compressible gas dynamics and MHD, the performance of the blending of more than one method is highly dependent on flow physics. Regardless of the choice of these blends of high-order methods, these test cases illustrate the improved stability and accuracy over the standalone single high-order method approach. Nine eighth-order structure-preserving methods were compared for a wide range of test cases in order to illustrate the effectiveness of the blending of the selected numerical method approaches. These current trends in the adaptive blending of high-order structure-preserving non-dissipative FDM with high-resolution high-order shock-capturing methods are also applicable to FVM, DG, FEM, and spectral element methods.

In addition, our studies have indicated that these current approaches can also improve stability and accuracy in rapidly developed unsteady flows containing steep gradients, as well as shear and shock waves, with the benefit of minimizing the added numerical dissipations.