# Local Characteristic Decomposition Based Central-Upwind Scheme

We propose novel less diffusive schemes for conservative one- and two-dimensional hyperbolic systems of nonlinear partial differential equations (PDEs). The main challenges in the development of accurate and robust numerical methods for the studied systems come from the complicated wave structures, such as shocks, rarefactions and contact discontinuities, arising even for smooth initial conditions. In order to reduce the diffusion in the original central-upwind schemes, we use a local characteristic decomposition procedure to develop a new class of central-upwind schemes. We apply the developed schemes to the one- and two-dimensional Euler equations of gas dynamics to illustrate the performance on a variety of examples. The obtained numerical results clearly demonstrate that the proposed new schemes outperform the original central-upwind schemes.

• 2 publications
• 1 publication
• 19 publications
• 3 publications
• 9 publications
03/05/2021

### Robust and accurate central algorithms for Multi-Component mixture equations with Stiffened gas EOS

Simple and robust algorithms are developed for compressible Euler equati...
03/24/2020

### Novel, simple and robust contact-discontinuity capturing schemes for high speed compressible flows

The nonlinear convection terms in the governing equations of compressibl...
08/07/2021

### Optimal parameters for numerical solvers of PDEs

In this paper we introduce a procedure for identifying optimal methods i...
05/10/2021

### Scalar auxiliary variable approach for conservative/dissipative partial differential equations with unbounded energy

In this paper, we present a novel investigation of the so-called SAV app...
04/08/2020

### Linearized Implicit Methods Based on a Single-Layer Neural Network: Application to Keller-Segel Models

This paper is concerned with numerical approximation of some two-dimensi...
05/04/2022

### Approximations of dispersive PDEs in the presence of low-regularity randomness

We introduce a new class of numerical schemes which allow for low regula...
11/25/2020

### A combination of Residual Distribution and the Active Flux formulations or a new class of schemes that can combine several writings of the same hyperbolic problem: application

We show how to combine in a natural way (i.e. without any test nor switc...

## 1 Introduction

This paper focuses on numerical solutions of hyperbolic systems of conservation laws, which in the two-dimensional (2-D) case, read as

 Ut+F(U)x+G(U)y=0, (1.1)

where and are spatial variables, is the time,

is a vector of unknowns,

and are the - and -flux functions, respectively.

It is well-known that the solutions of (1.1) may develop complicated wave structures including shocks, rarefactions, and contact discontinuities even for infinitely smooth initial data, and thus developing highly accurate and robust shock-capturing numerical methods for solving (1.1) is a challenging task.

Since the pioneering works [7, 22, 9], a large number of various methods had been introduced; see, e.g., the monographs and review papers [12, 23, 15, 38, 40, 2]

and references therein. Here, we focus on finite-volume (FV) methods, in which solutions are realized in terms of their cell averages and evolved in time according to the following algorithm. First, a piecewise polynomial interpolant is reconstructed out of the given cell averages, and then the evolution step is performed using the integral form of (

1.1). To this end, a proper set of time-space control volumes has to be selected. Depending on this selection, one may distinguish between two basic classes of finite-volume methods: upwind and central schemes. In upwind schemes, the spatial part of the control volumes coincides with the FV cells and therefore, one needs to (approximately) solve (generalized) Riemann problems at every cell interface; see, e.g., [9, 40, 2] and references therein. This helps the upwind schemes to achieve very high resolution. At the same time, it might be quite complicated and even impossible to solve the (generalized) Riemann problems for general systems of conservation laws. Central schemes offer an attractive simple alternative to the upwind ones. In staggered central schemes, first proposed for one-dimensional (1-D) systems in [30] and then extended to higher order [24, 29] and multiple number of dimensions [25, 13, 1], the control volumes use the staggered spatial parts so that the FV cell interfaces remain inside the control volumes. This allows to avoid even approximately solving any Riemann problems, which makes staggered central schemes easy to implement for a wide variety of hyperbolic systems. The major drawback of staggered central schemes is, however, their relatively large numerical dissipation as they basically average over the Riemann fans rather than resolving them.

In order to reduce the amount of excessive numerical dissipation present in central schemes, a class of central-upwind (CU) schemes have been proposed in [17, 19]. These schemes are based on nonuniform control volumes, whose spatial size taken to be proportional to the local speeds of propagation. This allows one to minimize the area over which the solution averaging occurs still without (approximately) solving any (generalized) Riemann problems. The upwind features of the CU schemes can be seen, for example, when they are applied to simpler systems. For instance, the CU scheme from [17]

reduces to the upwind one when it is applied to a system whose Jacobian contains only positive (only negative) eigenvalues. Another advantage of the CU schemes is related to the fact that unlike the staggered central schemes, they admit a particularly simple semi-discrete form. This observation is the basis of the modifications of the CU schemes we propose in this paper.

Even though the CU schemes from [17, 19] are quite accurate, efficient and robust tools for a wide variety of hyperbolic systems, higher resolution of the numerical solutions can be achieved by further reducing numerical dissipation. This can be done in a number of different ways, for example: (i) by introducing a more accurate evolution procedure, which leads to a “built-in” anti-diffusion term [18]; (ii) by implementing a numerical dissipation switch to control the amount of numerical dissipation present in the CU schemes [16]

; (iii) by obtaining more accurate estimates for the one-sided local speeds of propagation using the discrete Rankine-Hugoniot conditions

[8].

Another way to control the amount of numerical dissipation present in the CU schemes or any other FV methods is by adjusting the nonlinear limiting mechanism used in the piecewise polynomial reconstruction. It is well-known that sharper reconstructions may lead to larger numerical oscillations and a way to reduce these oscillations is to reconstruct the characteristic variables rather than the conservative ones; see, e.g., [33]. This can be done using the local characteristic decomposition (LCD); see, e.g., [3, 14, 28, 31, 33, 38, 42] and references therein.

In this paper, we modify the CU schemes from [17] by applying LCD to the numerical diffusion part of the schemes. The obtained new LCD-based CU schemes contain substantially smaller amount of numerical dissipation, which leads to a significantly improved resolution of the computed solution compared with the original CU schemes. As observed above, the key idea is that the new LCD-based CU scheme reduces to the upwind scheme when applied to 1-D linear hyperbolic systems

 Ut+AUx=0, (1.2)

where is a constant matrix (disregarding the sign of the eigenvalues of ). This feature suggests that the new CU schemes have more built-in upwinding compared with their predecessors.

The paper is organized as follows. In §2, we briefly describe the 1-D second-order FV CU scheme from [17]. In §3, we introduce the proposed new 1-D LCD-based CU scheme and show that the developed 1-D scheme reduces to the upwind scheme when applied to the linear system (1.2). In §4, we construct the 2-D LCD-based CU scheme. Finally, in §5, we test the proposed schemes on a number of 1-D and 2-D numerical examples for the Euler equations of gas dynamics. We demonstrate high accuracy, efficiency, stability, and robustness of the new LCD-based CU schemes, which clearly outperform the second-order CU scheme from [17, 20].

## 2 1-D Central-Upwind Scheme: A Brief Overview

In this section, we consider the 1-D hyperbolic system of conservation laws

 Ut+F(U)x=0, (2.1)

and describe the second-order semi-discrete CU scheme from [17]. To this end, we assume that the computational domain is covered with the uniform cells of size centered at and denote by cell averages of over the corresponding intervals , that is,

 \vbox{\hrule height 0.5pt width 100%\kern 1.72pt\hbox{\kern-0.5ptU\kern 0.0pt}}j(t):≈1Δx∫CjU(x,t)dx.

We also assume that at certain time , the cell average values are available and from here on we suppress the time-dependence of all of the indexed quantities for the sake of brevity.

According to [17], the computed cell averages (2.1

) are evolved in time by solving the following system of ordinary differential equations (ODEs):

 d\vbox{\hrule height 0.5pt width 100%\kern 1.72pt\hbox{\kern% -0.5ptU\kern 0.0pt}}jdt=−Fj+12−Fj−12Δx, (2.2)

where are the CU numerical fluxes given by

 Fj+12=a+j+12F−j+12−a−j+12F+j+12a+j+12−a−j+12+a+j+12a−j+12a+j+12−a−j+12(U+j+12−U−j+12). (2.3)

Here, and are the right/left-sided values reconstructed out of the given set of cell averages . The one-sided local speeds of propagation are estimated using the largest and the smallest eigenvalues of the Jacobian , which we denote by . This can be done, for example, by taking

 a+j+12=max{λd(A(U+j+12)),λd(A(U−j+12)),0}, (2.4) a−j+12=min{λ1(A(U+j+12)),λ1(A(U−j+12)),0}.

The (formal) order of the scheme (2.2)–(2.4) is determined by the accuracy the point values are reconstructed with and the order of the ODE solver used to integrate the ODE system (2.2) in time. In this paper, we use the LCD-based second-order piecewise linear reconstruction described in §2.1 and the three-stage third-order strong stability preserving (SSP) Runge-Kutta solver; see, e.g., [10, 11].

### 2.1 LCD-Based Piecewise Linear Reconstruction

In order to ensure at least second-order of accuracy of the CU scheme (2.2)–(2.4), one has to use a second-order piecewise linear reconstruction, which will be non-oscillatory provided the numerical derivatives are computed using a nonlinear limiter. A library of different limiters are available; see, e.g., [26, 30, 39, 2, 12, 23, 40]

. The limiters can be classified as dissipative, compressive or overcompressive; see

[26]. When applied to the conservative variables , dissipative limiters may introduce excessive numerical dissipation, which may results in oversmeared solution discontinuities. On the other hand, compressive and overcompressive limiters tend to lead to quite large oscillations near discontinuities and artificial sharpening smooth parts of the solution (compressive limiters may lead to kinks, that is, jumps in the first derivatives, while the overcompressive limiters may lead to jumps in the solution itself); see, e.g., [26]. Part of these difficulties can be overcome by reconstructing local characteristic variables as it was demonstrated in [33] in the context of higher-order WENO reconstructions.

In this paper, we use the minmod limiter, which is, according to [26], the most compressive out of dissipative limiters. As this limiter typically leads to quite large oscillations when applied to the conservative variables, we implement it in the LCD framework. Specifically, we first introduce the matrix , where is the Jacobian and is either a simple average or another type of average of the and states.

As long as the system (1.1) is strictly hyperbolic, we compute the matrices and such that is a diagonal matrix and introduce the local characteristic variables in the neighborhood of :

 Γk=R−1j+12\vbox{\hrule height 0.5pt width 10% 0%\kern 1.72pt\hbox{\kern-0.5ptU\kern 0.0pt}}k,k=j−1,…,j+2.

Equipped with the values , , and , we reconstruct in the cell by computing

 (Γx)j=minmod(2Γj−Γj−1Δx,Γj+1−Γj−12Δx,2Γj+1−ΓjΔx), (2.5)

and

 (Γx)j+1=minmod(2Γj+1−ΓjΔx,Γj+2−Γj2Δx,2Γj+2−Γj+1Δx), (2.6)

where the minmod function, defined as

 minmod(z1,z2,…):=⎧⎪⎨⎪⎩minj{zj}if zj>0∀j,maxj{zj}if zj<0∀j,0otherwise, (2.7)

is applied in the component-wise manner. The slopes (2.5) and (2.6) allow to evaluate

 Γ−j+12=Γj+Δx2(Γx)jandΓ+j+12=Γj+1−Δx2(Γx)j+1,

and obtain the corresponding point values of by

 U±j+12=Rj+12Γ±j+12. (2.8)
###### Remark 2.1

In Appendix A, we provide a detailed explanation on how the average matrix and the corresponding matrices and are computed in the case of the Euler equation of gas dynamics.

## 3 1-D LCD-Based Central-Upwind Scheme

In this section, we introduce a new 1-D LCD-based CU scheme, in which the amount of numerical dissipation is substantially reduced compared with the original CU scheme (2.2)–(2.3). To this end, we first rewrite the numerical flux (2.3) of the original CU scheme in the following form:

 Fj+12=Fj+Fj+12+Dj+12, (3.1)

where and is the numerical diffusion given by

 Dj+12 =a+j+12a+j+12−a−j+12[F−j+12−Fj+Fj+12]−a−j+12a+j+12−a−j+12[F+j+12−Fj+Fj+12] +a+j+12a−j+12a+j+12−a−j+12(U+j+12−U−j+12),

which, in turn, can be rewritten with the help of the matrix introduced in §2.1 and using (2.8) as follows:

 Dj+12 =Rj+12Pj+12R−1j+12[F−j+12−Fj+Fj+12]+Rj+12Mj+12R−1j+12[F+j+12−Fj+Fj+12] (3.2) +Rj+12Qj+12(Γ+j+12−Γ−j+12).

Here, , , and are the diagonal matrices

 (3.3) Qj+12=diag((Q1)j+12,…,(Qd)j+12)

with

 ((Pi)j+12,(Mi)j+12,(Qi)j+12)=1a+j+12−a−j+12(a+j+12,−a−j+12,a+j+12a−j+12),i=1,…,d. (3.4)

The main idea of the new the 1-D LCD-based CU scheme is to replace the constant entries (3.4) in the diagonal matrices , , and with the corresponding characteristic entries of the diagonal matrix . This is motivated by the fact that in the linear case with , the local characteristic speed might be different for each (diagonalized) component . Hence, a diffusion that depends on the corresponding local speed of the each component may lead to a sharper resolution of possible discontinuities.

Proposed modifications lead to the semi-discrete scheme

 d\vbox{\hrule height 0.5pt width 100%\kern 1.72pt\hbox{\kern% -0.5ptU\kern 0.0pt}}jdt=−FLCDj+12−FLCDj−12Δx, (3.5)

where the numerical fluxes are the LCD modifications of (3.1):

 FLCDj+12=Fj+Fj+12+DLCDj+12, (3.6)

with the following LCD modification of the numerical diffusion term (3.2)–(3.4):

 DLCDj+12 =Rj+12PLCDj+12R−1j+12[F−j+12−Fj+Fj+12]+Rj+12MLCDj+12R−1j+12[F+j+12−Fj+Fj+12] (3.7) +Rj+12QLCDj+12(Γ+j+12−Γ−j+12).

Here,

 PLCDj+12=diag((PLCD1)j+12,…,(PLCDd)j+12),MLCDj+12=diag((MLCD1)j+12,…,(MLCDd)j+12), QLCDj+12=diag((QLCD1)j+12,…,(QLCDd)j+12)

with

 ((PLCDi)j+12,(MLCDi)j+12,(QLCDi)j+12) (3.8) =⎧⎪ ⎪⎨⎪ ⎪⎩1(λ+i)j+12−(λ−i)j+12((λ+i)j+12,−(λ−i)j+12,(λ+i)j+12(λ−i)j+12)if (λ+i)j+12−(λ−i)j+12>ε,0otherwise,

and

 (λ+i)j+12 (λ−i)j+12

for .

Finally, in (3.8) is a very small desingularization constant, taken in all of the numerical examples reported in §5.

It is crucial to note that the proposed 1-D LCD-based CU scheme reduces to the second-order semi-discrete upwind scheme when applied to a linear hyperbolic system (1.2) with constant coefficients. This is proven in the following lemma.

###### Lemma 3.1

If

 F(U)=AU, (3.9)

where is a constant matrix, then the scheme (3.5)–(3.7) becomes the second-order semi-discrete upwind scheme.

Proof.

We note that in the linear case with constant coefficients, the corresponding matrix composed of the right eigenvectors of

is also a constant matrix, that is, , and matrices and reduce to

 PLCDj+12≡P=diag(max{sign(λ1),0},…,max{sign(λd),0}), (3.10) MLCDj+12≡M=diag(min{sign(λ1),0},…,min{sign(λd),0}),

while the matrix as .

Taking into account the above simplifications and substituting (3.9), (3.10) into (3.7), yields the following expression for the numerical diffusion term:

 DLCDj+12=RPR−1A[U−j+12−\vbox{\hrule height 0.5pt\kern 1.72pt\hbox{\kern-0.5% ptU\kern 0.0pt}}j+\vbox{\hrule height 0.5pt\kern 1.72pt\hbox{% \kern-0.5ptU\kern 0.0pt}}j+12]−RMR−1A[U+j+12−\vbox{\hrule height 0.5pt\kern 1.72pt\hbox{\kern-0.% 5ptU\kern 0.0pt}}j+\vbox{\hrule height 0.5pt\kern 1.72pt\hbox% {\kern-0.5ptU\kern 0.0pt}}j+12] (3.11) =RPR−1AR[R−1(U−j+12−\vbox{\hrule height 0.5pt\kern 1.72pt\hbox{\kern-0.5ptU\kern 0.0% pt}}j+\vbox{\hrule height 0.5pt\kern 1.72pt\hbox{\kern-0.5ptU% \kern 0.0pt}}j+12)]−RMR−1AR[R−1(U+j+12−\vbox{\hrule height 0.5pt\kern 1.72pt\hbox{\kern-0% .5ptU\kern 0.0pt}}j+\vbox{\hrule height 0.5pt\kern 1.72pt% \hbox{\kern-0.5ptU\kern 0.0pt}}j+12)] =RΛ+(Γ−j+12−Γj+Γj+12)+RΛ−(Γ+j+12−Γj+Γj+12),

where

 Λ+ =PR−1AR=PΛ=diag(max{λ1,0},…,max{λd,0}), Λ− =MR−1AR=MΛ=diag(min{λ1,0},…,min{λd,0}).

Finally, combining (3.9) and (3.11) together and including them into (3.6), we obtain

 FLCDj+12 =12A(\vbox{\hrule height 0.5pt\kern 1.72% pt\hbox{\kern-0.5ptU\kern 0.0pt}}j+\vbox{\hrule height 0.5pt% \kern 1.72pt\hbox{\kern-0.5ptU\kern 0.0pt}}j+1)+RΛ+(Γ−j+12−Γj+Γj+12)+RΛ−(Γ+j+12−Γj+Γj+12) (3.12) =12(A++A−)(\vbox{% \hrule height 0.5pt\kern 1.72pt\hbox{\kern-0.5ptU\kern 0.0pt}}j+\vbox{\hrule height 0.5pt\kern 1.72pt\hbox{\kern-0.5ptU\kern 0.0% pt}}j+1)+A+(U−j+12−Uj+Uj+12)+A−(U+j+12−Uj+Uj+12) =A+U−j+12+A−U+j+12,

where . This confirms that the scheme (3.5), (3.12) is the second-order semi-discrete upwind scheme.

## 4 2-D LCD-Based Central-Upwind Scheme

In this section, we generalize the 1-D LCD-based CU scheme introduced in §3 for the 2-D hyperbolic system of conservation laws (1.1). We design the 2-D LCD-based CU scheme in a “dimension-by-dimension” manner, so that it reads as

 d\vbox{\hrule height 0.5pt width 100%\kern 1.72pt\hbox{\kern% -0.5ptU\kern 0.0pt}}j,kdt=−FLCDj+12,k−FLCDj−12,kΔx−GLCDj,k+12−GLCDj,k−12Δy,

where

 FLCDj+12,k=Fj,k+Fj+1,k2+DLCDj+12,k,GLCDj,k+12=Gj,k+Gj,k+12+DLCDj,k+12.

Here, , , and and are the numerical diffusion terms defined by

 DLCDj+12,k =Rj+12,kPLCDj+12,kR−1j+12,k[FEj,k−Fj,k+Fj+1,k2] +Rj+12,kMLCDj+12,kR−1j+12,k[FWj+1,k−Fj,k+Fj+1,k2]+Rj+12,kQLCDj+12,k(ΓWj+1,k−ΓEj,k), DLCDj,k+12 =Rj,k+12PLCDj,k+12R−1j,k+12[GNj,k−Gj,k+Gj,k+12] −Rj,k+12MLCDj,k+12R−1j,k+12[GSj,k+1−Gj,k+Gj,k+12]+Rj,k+12QLCDj,k+12[ΓSj,k+1−ΓNj,k].

The matrices , and , are the matrices such that and are diagonal. Here, , with , , and , are either simple averages , or another type of averages of , and , states, respectively. The numerical fluxes and are defined by , , and the details of reconstructing the point values and , and for the 2-D Euler equations of gas dynamics are provided in Appendix B. Finally, the diagonal matrices , , , , , and are defined by

 PLCDj+12,k=diag((PLCD1)j+12,k,…,(PLCDd)j+12,k), PLCDj,k+12=diag((PLCD1)j,k+12,…,(PLCDd)j,k+12), MLCDj+12,k=diag((MLCD1)j+12,k,…,(MLCDd)j+12,k), MLCDj,k+12=diag((MLCD1)j,k+12,…,(MLCDd)j,k+12), QLCDj+12,k=diag((QLCD1)j+12,k,…,(QLCDd)j+12,k), QLCDj,k+12=diag((QLCD1)j,k+12,…,(QLCDd)j,k+12),

where for all

 ((PLCDi)j+12,k,(MLCDi)j+12,k,(QLCDi)j+12,k) =⎧⎪ ⎪⎨⎪ ⎪⎩1Δ(λi)j+12,k((λ+i)j+12,k,−(λ−i)j+12,k,(λ+i)j+12,k(λ−i)j+12,k)if Δ(λi)j+12,k>ε,0otherwise, ((PLCDi)j,k+12,(MLCDi)j,k+12,(QLCDi)j,k+12) =⎧⎪ ⎪⎨⎪ ⎪⎩1Δ(μi)j,k+12((μ+i)j,k+12,−(μ