A space-time smooth artificial viscosity method for nonlinear conservation laws

# A space-time smooth artificial viscosity method for nonlinear conservation laws

J. Reisner Los Alamos National Lab, XCP-4 MSF605, Los Alamos, NM 87544 J. Serencsa Department of Mathematics, UC Davis, One Shields Ave., Davis, CA 95616 S. Shkoller Department of Mathematics, UC Davis, One Shields Ave., Davis, CA 95616
###### Abstract

We introduce a new methodology for adding localized, space-time smooth, artificial viscosity to nonlinear systems of conservation laws which propagate shock waves, rarefactions, and contact discontinuities, which we call the -method. We shall focus our attention on the compressible Euler equations in one space dimension. The novel feature of our approach involves the coupling of a linear scalar reaction-diffusion equation to our system of conservation laws, whose solution is the coefficient to an additional (and artificial) term added to the flux, which determines the location, localization, and strength of the artificial viscosity. Near shock discontinuities, is large and localized, and transitions smoothly in space-time to zero away from discontinuities. Our approach is a provably convergent, spacetime-regularized variant of the original idea of Richtmeyer and Von Neumann, and is provided at the level of the PDE, thus allowing a host of numerical discretization schemes to be employed.

We demonstrate the effectiveness of the -method with three different numerical implementations and apply these to a collection of classical problems: the Sod shock-tube, the Osher-Shu shock-tube, the Woodward-Colella blast wave and the Leblanc shock-tube. First, we use a classical continuous finite-element implementation using second-order discretization in both space and time , FEM-C. Second, we use a simplified WENO scheme within our -method framework, WENO-C. Third, we use WENO with the Lax-Friedrichs flux together with the -equation, and call this WENO-LF-C. All three schemes yield higher-order discretization strategies, which provide sharp shock resolution with minimal overshoot and noise, and compare well with higher-order WENO schemes that employ approximate Riemann solvers, outperforming them for the difficult Leblanc shock tube experiment.

## 1 Introduction

### 1.1 Smoothing conservation laws

The initial-value problem for a general nonlinear system of conservation laws can be written as an evolution equation,

 ∂tU(x,t)+divF(U(x,t))=0 with U|t=0=U0, (1)

for an -vector defined on ()-dimensional space-time. Such partial differential equations (PDE) are both ubiquitous and fundamental in science and engineering, and include the compressible Euler equations of gas dynamics, the magneto-hydrodynamic (MHD) equations modeling ionized plasma, the elasticity equations of solid mechanics, and numerous related physical systems which possess complicated nonlinear wave interactions.

It is well known that solutions of (1) can develop finite-time shocks, even when the initial data is smooth, in which case, discontinuities of are propagated according to the so-called Rankine-Hugoniot conditions (see Section 2.1 below). It is important to develop stable and robust numerical algorithms which can approximate weak shock-wave solutions. Even in one-space dimension, nonlinear wave interaction such as two shock waves colliding, is a difficult problem when considering accuracy, stability and monotonicity. The challenge is maintaining higher-order accuracy away from the shock while approximating the discontinuity in an order- smooth transition region where denotes the spatial grid size.

As we describe below, a variety of clever discretization schemes have been developed and employed, particularly in one-space dimension, to approximate discontinuous solution profiles in an essentially non-oscillatory (ENO) fashion. These include, but are not limited to, total variation diminishing (TVD) schemes, flux-corrected transport (FCT) schemes, weighted essentially non-oscillatory (WENO) schemes, discontinuous Galerkin methods, artificial diffusion methods, exact and approximate Riemann solvers, and a host of variants and combinations of these techniques.

We develop a robust parabolic-type regularization of (1), which we refer to as the -method, which couples a modified set of equations for with an additional linear scalar reaction-diffusion equation for a new scalar field . Thus, instead of (1) we consider a system of equations, which use the solution as a coefficient in a carefully chosen modification of the flux. As we describe in detail below, the solution is highly localized in regions of discontinuity, and transitions smoothly (in both and ) to zero in regions wherein the solution is smooth. Further, as , we recover the original hyperbolic nonlinear system of conservation laws (1).

### 1.2 Numerical discretization

In the case of 1-D gas dynamics, the construction of non-oscillatory, higher-order, numerical algorithms such as ENO by Harten, Engquist, Osher & Chakravarthy Harten1987231 () and Shu & Osher Shu1988439 (), Shu198932 (); WENO by Liu, Osher, & Chan LiOsCh1994 () and Jiang & Shu Jiang1996202 (); MUSCL by Van Leer VanLeer1979101 (), Colella Colella1985104 (), and Huynh Huynh19951565 (); or PPM by Colella & Woodward Colella1984174 () requires carefully chosen reconstruction and numerical flux.

Such numerical methods evolve cell-averaged quantities; to calculate an accurate approximation of the flux at cell-interfaces, these schemes reconstruct th-order () polynomial approximations of the solution (and hence the flux) from the computed cell-averages, and thus provide th-order accuracy away from discontinuities. See, for example, the convergence plots of Greenough & Rider Greenough2004259 () and Liska & Wendroff Liska2003995 (). Given a polynomial representation of the solution, a strategy is chosen to compute the most accurate cell-interface flux, and this is achieved by a variety of algorithms. Centered numerical fluxes, such as Lax-Friedrichs, add dissipation as a mechanism to preserve stability and monotonicity. On the other hand, characteristic-type upwinding based upon exact (Godunov) or approximate (Roe, Osher, HLL, HLLC) Riemann solvers, which preserve monotonicity without adding too much dissipation, tend to be rather complex and PDE-specific; moreover, for strong shocks, other techniques may be required to dampen post-shock oscillations or to yield entropy-satisfying approximations (see Quirk Quirk1994555 ()). Again, we refer the reader to the papers Greenough2004259 (), Liska2003995 () or Colella & Woodward Colella1984115 () for a thorough overview, as well as a comparison of the effectiveness of a variety of competitive schemes.

Majda & Osher Majda1977671 () have shown that any numerical scheme is at best, first-order accurate in the presence of shocks or discontinuities. The use of higher-order numerical schemes is, nevertheless, imperative for the elimination of error-terms in the Taylor expansion (in mesh-size) and the subsequent limiting of truncation error. Moreover, higher-order schemes tend to be less dissipative than there lower-order counterparts, as discussed by Greenough & Rider Greenough2004259 (); therein, a comparison between a nd-order PLMDE scheme and a th-order WENO scheme demonstrates the improved resolution of intricate fine structure afforded by th-order WENO, while simultaneously providing far less clipping of local extrema than PLMDE.

In multi-D, similar tools are required to obtain non-oscillatory numerical schemes, but the multi-dimensional analogues to those described above are generally limited by mesh considerations. For structured grids (such as products of uniform 1-D grids), dimensional splitting is commonly used, decomposing the problem into a sequence of 1-D problems. This technique is quite successful, but stringent mesh requirements prohibits its use on complex domains. Moreover, applications to PDE outside of variants of the Euler equations may be somewhat limited. For further discussion of the limitations of dimensional splitting, we refer the reader to Crandall & Majda Crandall80 (), and Jiang & Tadmor Jiang98 (). For unstructured grids, dimensional splitting is not available and alternative approaches must be employed, necessitated by the lack of multi-D Riemann solvers. WENO schemes on unstructured triangular grids have been developed in Hu & Shu Hu199997 (), but using simplified methods, which employ reduced characteristic decompositions, can lead to a loss of monotonicity and stability.

Algorithms that explicitly introduce diffusion provide a simple way to stabilize higher-order numerical schemes and subsequently remove non-physical oscillations near shocks. In the mathematical analysis of conservation laws (and in the truncation error of certain discretization schemes), the simplest parabolic-regularization is by the addition of a uniform linear viscosity. Choosing a constant , which depends upon mesh-size and sometimes velocity or wave-speed, and adding

 β(Δx)∂2xU(x,t) (2)

to the right hand side of (1) provides a uniformly parabolic regularization of the hyperbolic conservation laws, and its discrete implementation smears sharp discontinuities across -regions and thus adds stabilization, but unfortunately, at the cost of accuracy. With the addition of uniform linear viscosity, shocks and discontinuities are captured in a non-oscillatory fashion, but the transition region form left to right state, which approximates the discontinuity, tends to grow over time. Moreover, since viscosity is applied uniformly over the entire domain , the benefits of a higher-order scheme (away from the discontinuity) may be lost, and the accuracy may often reduce to merely first-order. In practical implantation in a numerical scheme, the use of viscosity should be localized in regions of shock (so as to stabilize the scheme), limited at contact discontinuities (to avoid over-smearing the sharp transition), and very small in smooth regions away from discontinuities. Achieving these requirements allows higher-order approximation of smooth flow and sharp, non-oscillatory, resolution of shocks and discontinuities. Naturally, this necessitates that the amount of added viscosity be a function of the solution.

The pioneering papers of Richtmyer Richtmyer1948 (), Von Neumann & Richtmyer VonNeumann1950232 (), Lax & Wendroff Lax1960217 (), and Lapidus Lapidus1967154 () suggest the introduction of nonlinear artificial viscosity to equations (1) in a form similar to the following expression:

 β(Δx)2∂x(|∂xu(x,t)|∂xU(x,t)). (3)

We refer the reader to the classical papers of Gentry, Martin, & Daly GeMaDa1966 () and Harlow & Amsden HaAm1971 () for an interesting discussion on artificial viscosity. Specifically, Gentry, Martin, & Daly GeMaDa1966 () define the nonlinear viscosity of the type (3) to be artificial viscosity, and show that the linear viscosity (2), scaled by the magnitude of local velocity, arises as truncation error (in finite-difference approximations). The latter is responsible for stabilizing the transport of sound waves, while (3) stabilizes the steepening of sound waves.111We are indebted to the anonymous referee for clarifying this point for us.

We are primarily concerned with the steepening of sound waves, and shall term artificial viscosity of the type (3) as classical artificial viscosity. Formally, the use of (3) produces the required amount of viscosity near shocks but allows for second-order accuracy in smooth regions. On the other hand, the diffusion coefficient is precisely the quantity which loses regularity (or smoothness) near shock discontinuities. Also, the constant must be larger than one to control numerical oscillations behind the shock wave, which in turn overly diffuses the waves and produces incorrect wave speeds.

Alternative procedures have been proposed. For streamline upwind Petrov-Galerkin schemes (SUPG), Hughes & Mallet Hughes1986329 () and Shakib, Hughes, & Johan Shakib1991141 () use residual-based artificial viscosity. Guermond & Pasquetti Guermond2008801 () present a similar, entropy-residual-based scheme for use in spectral methods. Persson & Peraire Persson2006112 () develop a method based upon decay of local interpolating polynomials for discontinuous Galerkin schemes. Later, Barter & Darmofal Barter20101810 () use a reaction-diffusion equation to provide a regularized variant of this approach.

Our approach is similar to Barter20101810 () in that it uses a reaction-diffusion equation to calculate a smooth distribution of artificial viscosity. Instead of regularizing a DG-based noise-indicator that allows for the growth of viscosity near shocks, we regularize the classical artificial viscosity of Lapidus1967154 (), using a gradient based approach for this source term. This approach yields both a discretization- and PDE-independent methodology which can be generalized to multiple dimensions by regularizing a similar viscosity to that in Löhner, Morgan, & Peraire Lohner1985141 ().

In 1-D, our approach proves to be a simple way of circumventing the need for characteristic or other a priori information of the exact solution to remove oscillations in higher-order schemes. Due to the simple and discretization-independent nature of our method, we expect our methodology to be useful for a wide range of applications.

### 1.3 Outline of the paper

In Section 2, we introduce the -method for the compressible Euler equations in one space dimension. We show that the -method is Galilean invariant and that solutions of the -method converge to the unique weak solutions of the Euler equations in the limit of zero mesh size. We also show the relative smoothness of our new viscosity coefficient with respect to the classical artificial viscosity of Richtmyer and Von Neumann, and we demonstrate the ability of the -method to remove downstream oscillation in slowly moving shocks.

In Section 3, we give a brief outline of the numerical schemes whose solutions are used in this paper. First, we outline a second-order, continuous Galerkin finite-element method. Second, we outline a simple WENO-based finite-volume scheme, only upwinding via the sign of the velocity (no Riemann-solvers or characteristic decompositions in primitive variables). The resulting schemes applied to the -method are referred to as FEM-C and WENO-C, respectively. Third, we outline the central-finite-difference scheme of Nessyahu and Tadmor (NT), a simple scheme, easily generalizable to multi-D Nessyahu1990408 (). Like our FEM-C scheme, the NT-scheme is at best, second-order, and does not require specialized techniques for upwinding. Fourth, we outline a Godunov-type characteristic decomposition-based WENO scheme (WENO-G) developed by Rider, Greenough & Kamm rideretal () which utilizes a variant of a Godunov/Riemann-solver as upwinding, providing a very competitive scheme for modeling the collision of very strong shocks.

In Section 4, we consider the classical shock-tube problem of Sod. With the Sod shock problem, we apply our FEM-C scheme and compare with the classical viscosity approach. We then compare the FEM-C scheme with the two standalone methods, NT and WENO-G.

In Section 5, we consider the moderately difficult problem of Osher-Shu, modeling the interaction of a mild shock with an entropy wave. We compare FEM-C to NT and WENO-G in which the differences are more significant than in the Sod-shock comparisons. We show that WENO-C compares well with WENO-G; on the other hand, the simple WENO scheme without the -method and without the Gudonov-based characteristic solver also does well in modeling the Osher-Shu test case.

In Section 6, we consider the numerically challenging Woodward-Colella blast wave simulation, which models the collision of two strong interacting shock fronts. Though the FEM-C scheme performs better than NT, both second-order schemes are somewhat out-performed by the higher-order WENO-G method (with characteristic solver). On the other hand, WENO-C compares well with WENO-G, having slightly less damped amplitudes with the same shock resolution.

Finally, in Section 7, we consider the Leblanc shock-tube, an extremely difficult test case consisting of a very strong shock. For this problem, devise two strategies to demonstrate the use of the -method. In the first strategy, we use our simplified WENO-C scheme with a right-hand side term for the energy equation that relies on a second -equation which smooths gradients of . We obtain an excellent approximation of the notoriously difficult contact discontinuity for internal energy, while maintaining an accurate shock speed; simultaneously, we avoid generating large overshoots at the contact discontinuity, which would indeed occur without the use of the -method. For our second strategy, we show that WENO with the Lax-Friedrichs flux can be significantly improved with the addition of the -method. We call this algorithm WENO-LF-C, and show that by using just one -equation (as we have for all of the other test cases), we can sharply resolve the contact discontinuity for the internal energy, with accurate wave speed, and without overshoots.

## 2 The C-method

We begin with a description of the 1-D compressible Euler equations, written as a 3x3 system of conservation laws. We then explain our parabolic regularization, which we call the -method.

### 2.1 Compressible Euler equations

The compressible Euler equations set on a 1-D space domain , and a time interval are written in vector-form as the following coupled system of nonlinear conservation laws:

 ∂tu(x,t)+∂xF(u(x,t)) =0, x∈I,t>0, (4a) u(x,0) =u0(x), x∈I,t=0, (4b)

where the -vector and flux function are defined, respectively, as

 u=⎛⎜⎝ρmE⎞⎟⎠andF(u)=⎛⎜ ⎜ ⎜⎝mm2ρ+pmρ(E+p)⎞⎟ ⎟ ⎟⎠,

and

 u0(x)=⎛⎜⎝ρ0(x)m0(x)E0(x)⎞⎟⎠

denotes the initial data for the problem. The variables , and denote the density, momentum, and energy density of a compressible gas, while denotes the pressure function. It is necessary to choose an equation-of-state , and we use the ideal gas law, for which

 p=(γ−1)(E−m22ρ), (5)

where denotes the adiabatic constant. The equations (4) are indeed conservation laws, as they represent the conservation of mass, momentum, and energy in the evolution of a compressible gas. The velocity field is obtained from momentum and density via the identity

 u=mρ.

Inverting the relation (5), we see that the energy density is a sum of kinetic and potential energy density functions:

 E=ρu22kinetic+pγ−1potential.

The gradient (or Jacobian) of the flux vector is given by

 DF(u)=⎡⎢ ⎢ ⎢ ⎢⎣010(γ−3)m22ρ2(3−γ)mργ−1−γEmρ2+(γ−1)m3ρ3γEρ+(1−γ)3m22ρ2γmρ⎤⎥ ⎥ ⎥ ⎥⎦

with eigenvalues

 λ1=u+c, λ2=u, λ3=u−c, (6a)

where denotes the sound speed (see, for example, Toro Toro2009 ()). These eigenvalues determine the wave speeds.

The behavior of the various wave patterns is greatly influenced by the speed of propagation; as such, we define the maximum wave speed to be

 [S(u)](t)=maxi=1,2,3maxx∈I{|λi(x,t)|}. (6b)

as a function of .

We are interested in solutions with shock waves and contact discontinuities. The Rankine-Hugoniot (R-H) conditions determine the speed of the moving shock discontinuity, as well as the speed of a contact discontinuity. For a shock wave discontinuity, the R-H condition can be stated

 F(ul)−F(ur)=s(ul−ur)

where the subscript denotes the state to the left of the discontinuity, and the subscript denotes the state to the right of the discontinuity. In general, the following three jump conditions must hold:

 ml−mr =s(ρl−ρr) ((3−γ)m2l2ρ2l+(γ−1)El)−((3−γ)m2r2ρ2r+(γ−1)Er) =s(ml−mr) (γElmlρl−γ−12m3lρ2l) =s(El−Er).

There can be non-uniqueness for weak solutions that have jump discontinuities, unless entropy conditions are satisfied (see the discussion in Section 2.9.4). So-called viscosity solutions are known to satisfy the entropy condition (and are hence unique) and are defined as the limit as of a sequence of solutions to the following parabolic equation:

 ∂tuϵ+∂xF(uϵ) =ϵ∂xxuϵ, t>0, (7a) uϵ =u0, t=0. (7b)

In the isentropic setting, for bounded initial data with bounded variation, solutions converge to the unique entropy solution of (4) as (see DiPerna DiPerna1983 () and Lions, Perthame, & Souganidis LiPeSo1996 ()). For non-isentropic dynamics, the same result holds if the initial data has small total variation (see Bianchini & Bressan BiBr2005 ()). Moreover, if the initial data is regularized, then solutions to (7) are smooth in both space and time, and the discontinuity is approximated by a smooth function, transitioning from the left-state to the right-state over an interval whose length is .

Some of the classical finite-differencing schemes, such as the Lax-Friedrichs discretization, is dissipative to second-order and effectively behaves as a discrete version of (7). The uniform nature of such diffusion does not distinguish between discontinuous and smooth flow regimes, and thus adds unnecessary dissipation in regions of the wave profile which do not require any numerical stabilization. Such uniform dissipation contributes to a non-phyiscal damping of entropy waves, over-diffusion and smearing of contact discontinuities, and changes the wave speeds. Ultimately, uniform artificial viscosity is not ideal; rather, artificial viscosity should be added in a localized and smooth manner.

### 2.2 Classical artificial viscosity

The idea of adding localized artificial viscosity to capture discontinuous solution profiles in numerical simulations dates back to Richtmyer Richtmyer1948 (), Von Neumann & Richtmyer VonNeumann1950232 (), Lax & Wendroff Lax1960217 (), Lapidus Lapidus1967154 () and a host of other reseachers. The idea behind classical artificial viscosity is to refine the uniform viscosity on the right-hand side of equation (7a) with

 ∂tuϵ+∂xF(uϵ)=βϵ2∂x(|∂xuϵ|∂xuϵ),t>0, (8)

for a suitably chosen constant , which may depend upon the numerical discretization scheme.

When the velocity exhibits a jump discontinuity (i.e., at a shock), the quantity is ; however, away from shocks, where the velocity is smooth, remains uniformly bounded in , and in such smooth regions, (8) adds significantly less viscosity than (7a). On the other hand, as we shall demonstrate in Figure 1, the use of as a coefficient in the smoothing operator, can lead to spurious oscillations in the solution, caused by the lack of regularity in the quantity .

Formally, the use of the localizing coefficient corrects for the over-dissipation of the uniform viscosity in (7), and a variety of numerical methods have employed some variant of this idea, achieving methods that are nominally non-oscillatory near shocks while maintaining second-order accuracy away from shocks. However, as we have already noted, the quantity may become highly irregular near shock discontinuities, and may thus require setting the constant in order to stabilize incipient numerical oscillations (see Section 4 for evidence to this observation). While this increase in does not effect the asymptotic accuracy of the scheme, it is clearly beneficial to take as small as possible to preserve the correct amplitude and wave speed.

The loss of regularity of the coefficient suggests that a smoothed version of would greatly benefit the dynamics. Smoothing in space is not sufficient, as we must ensure smoothness in time as well. Hence, we propose our -method, which indeed provides a regularized version of (8) and allows for the use of much smaller values of (less localized artificial dissipation), higher accuracy, and practical viability.

### 2.3 C-method for compressible Euler

Analogous to (8), we control the amount of viscosity in (7a) by the use of a function of space and time. This function is a solution to a reaction-diffusion equation, coupling to the evolution of . The mechanism of diffusion, smoothing/spreading the sharp peaks localized around discontinuities in the velocity competes with the mechanism of reaction, pushing to zero. Subsequently, our -method yields a smooth, yet sharp, distribution of artificial viscosity yielding regularized similar to (7) on which we can build high-resolution numerical schemes.

For fixed we choose . Then, for each , we find

 uϵ(x,t)=⎛⎜⎝ρϵ(x,t)mϵ(x,t)Eϵ(x,t)⎞⎟⎠   and   Cϵ(x,t)

as solutions of the following parabolic system of (viscous) conservation laws:

 ∂tuϵ+∂xF(uϵ) =∂x(~βϵ2Cϵ,δ∂xuϵ), t>0, (9a) ∂tCϵ−S(uϵ)∂2xCϵ+S(uϵ)ϵCϵ =S(uϵ)G(∂xuϵ), t>0, (9b) (uϵ,Cϵ) =(uϵ0,G(∂xuϵ0)), t=0, (9c)

where for a fixed positive constant , and . The forcing to equation (9b) is defined as

 G(∂xuϵ)=|∂xuϵ|maxI |∂xuϵ| (10)

is defined by (6), and denotes a regularization of the initial data which we discuss below. We also note that the scaling factor in , given by , is included only to make comparisons with the classical artificial viscosity approach, but is in no way necessary.

### 2.4 Regularization of initial data for use with FEM-C

Unlike numerical algorithms which advance cell-averaged quantities, the finite-element method relies upon polynomial interpolation of nodal values, and requires solutions to be continuous across element boundaries in order for the interpolation to converge. As such, the use of discontinuous initial data produces Gibbs-type oscillations, at least on very short time intervals. To avoid this spurious behavior, it is advantageous to smooth discontinuous initial profiles.

More specifically, we provide a hyperbolic-tangent smoothing for initial data for our FEM-C scheme. Since pointwise evaluation is well-defined for smooth functions, the finite-element discretization scheme can interpolate the regularized data and generate appropriate initial states.

For an interval , we denote the indicator function

 1[a,b](x)={1,x∈[a,b],0,x∉[a,b], (11)

and consider initial conditions with components of the form

 (uϵ0(x))i=Li∑j=11[aij,bij](x)fij(x),

where the are pairwise-disjoint (in ) and each are smooth.

We then define the regularized initial condition

 (uϵ0(x))i=Li∑j=11ϵ[aij,bij](x)fij(x),

where

 1ϵIij(x)=12[tanh(x+bijϵ)−tanh(x+aijϵ)].

This regularization procedure achieves approximations of exponential-order away from discontinuities; near discontinuities, it is a first-order approximation, when measured in the -norm. Specifically, if is smooth in , then the -norm of the error

 ∥(u0)i−(uϵ0)i∥L1(ω)=∫ω∣∣(u0(x))i−(uϵ0(x))i∣∣ dx=O(ϵp) (12)

for any positive integer . Alternatively, if is discontinuous somewhere in , the -norm of the error

 ∥(ui0)−(uϵ0)i∥L1(Ω)=O(ϵ). (13)

These observations assert that our regularization of the initial data allows for higher-order approximation of the initial data and is analogous to the averaging procedure required by Majda & Osher Majda1977671 ().

### 2.5 A compressive modification of the forcing G in the C-equation

The function in (10) is chosen in such a manner so that is large where there are sharp transitions in the velocity field . In compressive regions (i.e., where or ), sharp transitions over lengths of correspond to shocks and artificial viscosity is required so that remains smooth. In expansive regions, corresponding to rarefactions, artificial viscosity is not generally necessary.

These observations motivate the following alternative forcing function:

 Gcomp(∂xuϵ)=|∂xuϵ|Imax |∂xuϵ|1(−∞,0)(∂xuϵ) (14)

where the indicator function introduces viscosity only in regions of compression.

The ability to use such a switch is heavily dependent on the use of a space-time smoothing. Since the velocity in many numerical schemes may become oscillatory near shocks, such a switch can become discontinuous between adjacent cells/elements. However, the space-time nature of the -equation resolves this issue, providing a smooth artificial viscosity profile.

This modified function typically increases accuracy in Euler simulations, but can lead to a loss of stabilization. For our FEM-C approach, where the stabilizing effects of artificial viscosity are necessary to dampen noise, the use of is restricted to the problems of Sod and Osher-Shu, which contain only moderately strong shocks.

### 2.6 Moving to the discrete level

The use of the C-equation yields smooth solutions and thus we expect that a variety of higher-order discretization techniques, with sufficiently small and , could provide accurate, non-oscillatory approximations. In our implementation, artificial viscosity spreads discontinuities over regions of size . Thus, given a particular initial condition, final time, discretization scheme, etc., we choose such that the scaling produces non-oscillatory profiles.

We also note that the initial condition for , given in (9c) is chosen so to guarantee the coefficients of diffusion in (9a) are smooth up to . Moreover, choosing such initial conditions for allows one to recover the classical artificial viscosity as . As stated, these initial conditions may require a smaller time-step (by a factor of ) for the first few time-steps. In practice, taking is an effective simplification to eliminate the need for smaller initial time-steps. Alternatively, we can solve an elliptic PDE for at the initial time and similarly eliminate that concern.

### 2.7 The C-method under a Galilean-transformation

We begin our discussion for the case of constant entropy. The Galilean invariance of the isentropic Euler equations results from the advective nature of the PDE. Since we solve a modified equation (coupled with the additional -equation) it is of interest to know to what extent Galilean invariance is preserved. For simplicity, we assume that

 p(x,t)=ρ(x,t)2.

(The choice corresponds to the shallow water equations, but any other choice of works in the same fashion.)

Given a fixed we define the change in independent variables

 ~x=x−vt,~t=t,

denoting and the analogous change in the dependent variables

 ~ρ(~x,~t)=ρ(~x+v~t,~t),~u(~x,~t)=u(~x+v~t,~t)−v. (15)

A simple calculation yields

 ∂~t~ρ+∂~x(~ρ~u) =[∂tρ+∂x(ρu)]∘ϕ, (16a) ∂~t(~ρ~u)+∂~x(~ρ~u2+~p) =[∂t(ρu)+∂x(ρu2)]∘ϕ+∂~x~p−v[∂tρ+∂x(ρu)]∘ϕ. (16b)

We further have that

 ~p(~x,~t)=p(~x+v~t,~t), (17)

so that the mass and momentum equations are, in fact, Galilean invariant in the absence of artificial viscosity.

With the -method employed, (16) transforms to

 ∂~t~ρ+∂~x(~ρ~u) =[∂x(C∂xρ)]∘ϕ, (18a) ∂~t(~ρ~u)+∂~x(~ρ~u2+~p) =(∂x{C∂x[ρ(u−v)]})∘ϕ, (18b)

where we let , and drop the superscript for notational convenience.

Examining (9b), we see that the equation for is not Galilean invariant, but this is not a physical quantity, but can rather be viewed as a parameter to the modified system of conservation laws. As such we define the behavior of under Galilean transformations as follows:

 ~C(~x,~t)=C(~x+v~t,~t).

With this definition of , we find that

 ∂~t~ρ+∂~x(~ρ~u) =[∂~x(~C∂~x~ρ)], (19a) ∂~t(~ρ~u)+∂~x(~ρ~u2+~p) ={∂~x[~C∂~x(~ρ~u)]}, (19b)

and hence the -method for isentropic Euler retains the Galilean invariance.

We remark that in the absence of artificial viscosity on the right-hand side of the mass equation, the artificial flux term in the momentum equation is modified according to (34) below. This modification ensure Galilean invariance when the mass equation is left unchanged, which is the strategy employed for our WENO-C scheme.

Next, since the Galilean symmetry is for the smooth solutions (for which classical derivatives are well-defined), and since smooth velocity fields simply transport the entropy function, it is thus a consequence of the transport of entropy, that Galilean invariance holds for the non-isentropic case as well. The important of a numerical approximation to capture the Galilean invariant solution is fundamental to the initiation of the Kelvin-Helmholtz instability and other basic instabilities present in the Euler equation; see Robertson, Kravtsov, Gnedin, and Rudd RoKrGnAbRu2010 () for a thorough discussion. In this connection, we next examine long wavelength instabilities which can arise for very slowly moving shock waves.

### 2.8 Regularization through the C-equation

It is of interest to examine the relative smoothness of to its unsmoothed counterpart , and to determine the effect of this smoothing relative to the classical artificial viscosity approach. In Figure 1 we provide two plots demonstrating the effect of the -method. In Figure 1(a) we see that the -equation provides a smoothened viscosity profile compared to the classical approach. Alternatively, in Figure 1(b) we plot using the compression-switch modification versus using purely the quantity (not smoothed by the -equation) as a viscosity. In both cases we see how the -method provides a far smoother profile with roughly the same magnitude as the non-smoothened approach.

A useful feature of the -method is the ability to tune parameters in the -equation to generate non-oscillatory behavior. Though we are quite explicit on the form of the -equation in (9b), a simple modification allows for the diffusion coefficient to be problem dependent, i.e. allowing for a choice of positive constant and replacing the diffusion term with

 −γS(uϵ)∂2xCϵ.

In most of the forthcoming experiments, we fix , but we note that choosing larger can yield smoother solution profiles as the profile of will be less localized. The parameter is a time-relaxation parameter, and can be viewed in an analogous fashion to the time-relaxation parameter present in Cahn-Hilliard and Ginzburg-Landau theories. For very slow moving shocks, the time-relaxation can be adjusted to scale with the shock speed.222We note that is inversely proportional to the Mach number and its precise functional relation shall be examined in future work.

We find this to be an effective approach for the flattening procedure discussed in Colella1984174 () for removing oscillations that form to the left of a slowly right-moving shock. Moreover, Roberts Roberts1990 () concludes that a differentiable form of the numerical flux construction appears necessary to remove downstream long-wavelength oscillations caused by slow shock motion. We use the -method to analyze this.

Using the slow-shock initial conditions outlined in Quirk Quirk1994555 (), in Figure 2 we show the success of the FEM-C (outlined below in Section 3.2) in removing these oscillations when choosing (Fig. 2(a)) and (Fig. 2(b)).

### 2.9 Convergence of the C-method in the limit of zero mesh size

#### 2.9.1 The isentropic case

We sketch the proof for the isentropic Euler equations given by

 (ρu)t+(ρu2+p)x =0, (20a) ρt+(ρu)x =0, (20b) p(ρ) =ργ, (20c)

where .

To simplify the notation, we set , and set the momentum . Following (9), we write the -method version of (20) as

 mϵt+[(mϵ)2/ρϵ+pϵ]x =ϵ2[Cϵ,δmϵx]x, (21a) ρϵt+mϵx =ϵ2(Cϵ,δρx)x, (21b) pϵ(ρϵ) =(ρϵ)γ, (21c) Cϵt−S(uϵ)Cϵxx+S(uϵ)ϵCϵ =S(uϵ)G(uϵx), (21d)

or (21a,b) can equivalently in terms of the vector and flux as

 uϵt+f(uϵ)x=ϵ2[Cuϵx]x, (21’)

where denotes a diagonal 2x2 matrix with entries which is strictly positive-definite. Recall that , satisfies , and that , with denoting the sound speed. On any time interval , the maximum wave speed is uniformly strictly positive; thus, as the initial data for , the maximum principle shows that must be non-negative. We remark that while the use of as the coefficient is not required for the numerics, as is taken much smaller than the mesh size , strict positivity of simplifies the proof of regularity of solutions to (21) as well as the convergence argument.

To avoid issues with spatial boundaries, we shall assume periodic boundary conditions for our spatial domain. Note that in this case, the fundamental theorem of calculus shows that and that mass is conserved.

#### 2.9.2 The basic energy law

In order to prove that solutions to (21) converge to solutions of (20), we must establish -independent estimates for solutions of (21). To do so, we multiply equation (21a) by , integrate over our spatial domain, and make use of the equation (21b) to find that any weak solution to (21) must verify the basic energy law

 ddt[∫12ρϵ(uϵ)2dx+1γ−1∫pϵdx]≤−ϵ2∫Cϵ,δρϵ(uϵx)2dx−ϵ2γ∫Cϵ,δ(ρϵ)γ−2(ρϵx)2dx. (22)

(The inequality in (22) is due to the lower semi-continuity of weak convergence and is replaced with equality for solutions which are sufficiently regular.) Thus, the total energy of isentropic gas dynamics is dissipated according to the right-hand side of (22), and for each , we see that and are square-integrable (in ) for almost every instant of time, if the density , that is, if avoids vacuum. We shall explain below that this is indeed the case.

#### 2.9.3 Regularity of solutions uϵ

The reaction-diffusion equation (21d) is a uniformly parabolic equation. According to the energy law, and as a consequence of Sobolev’s theorem, is a bounded function; furthermore, the right-hand side of (21d) is in . It is standard, from the regularity theory of uniformly parabolic equations, that then has two spatial (weak) derivatives which are square-integrable. This, in turn, shows that for , solutions possess three spatial (weak) derivative which are square-integrable for almost every instant of time. This implies that solutions are classically differentiable in both space and time.

Furthermore, by using the symmetrizing matrix we can show that are, independently of and , uniformly bounded in the Sobolev space (consisting of measurable functions with two weak derivatives in ), and thus we may take a pointwise limit of this sequence as , in the event that the time-interval is sufficiently small as to ensure that a shock has not yet formed. Of course, we are interested, in convergence to discontinuous profiles, so we address this next.

#### 2.9.4 Convergence to the entropy solution

We shall now provide a sketch of the limit as . A function is called an entropy for (20) with entropy flux if smooth solutions satisfy the additional conservation law

 η(u)t+q(u)x=0. (23)

In non-conservative form, (20) and (23) are written as

 ut+∇f(u)ux=0,  ∇η(u)ut+∇q(u)ux=0,

from which we obtain the compatibility condition between and ,

 ∇η(u)∇f(u)=∇q(u). (24)

The pair satisfy (23) if and only if condition (24) holds. Moreover, a weak solution to (20) is the unique entropy solution if

 η(u)t+q(u)x≤0. (25)

For isentropic gas dynamics we can set

 η(m,ρ)=m22ρ+ργγ−1

which is the total energy, with corresponding entropy flux

We observe that is strongly convex as long as .

For the sequence of solution of (21), suppose that as , converges boundedly (almost everywhere) to a weak solution of (20). We claim that if satisfy (23), then (25) holds in the distributional sense. To see that this is the case, we take the inner-product of with equation (21’), and find that

 η(uϵ)t+q(uϵ)x =ϵ2∇η(uϵ)[Cuϵx]x =ϵ2[Cη(uϵ)x]x−ϵ2[uϵx]TC∇2η(uϵ)uϵx.

Integrating over the spatial domain and then over the time interval yields

 ∫η(uϵ(x,T))dx−∫η(uϵ(x,0))dx=−ϵ2∫T0∫[uϵx]TC∇2η(uϵ)uϵxdxdt,

from which it follows that

 ∫T0∫|ϵuϵx|2dxdt≤¯c (26)

where the constant depends upon , the minimum value of density, and the entropy in the initial data. For a smooth, non-negative test function with compact support in the strip ,

 ∬η(uϵϕt+q(uϵ)ϕxdxdt=ϵ∬C(ϵuϵ)xϕxdxdt+∬ϵ2[uϵx]TC∇2η(uϵ)uϵxϕdxdt.

Thanks to (26), the first term on the right-hand side goes to zero like , while the second term is non-negative, since is positive-definite (since is strongly convex) as is . Thus, as , we recover the entropy inequality (25).

It remains to discuss the assumptions concerning the bounded convergence of to , as well as the uniform bound from below . The argument relies on finding a priori bounds on the amplitudes of solutions to (21). If it is the case that uniformly in ,

 |uϵ|≤M and 0<ν≤ρϵ,

then the compensated-compactness approach for isentropic Euler pioneered by DiPerna DiPerna1983 () and made much more general by Lions, Perthame, & Souganidis LiPeSo1996 () provides a subsequence of converging pointwise (almost everywhere) to a solution of (20).

For isentropic gas dynamics, our approximation (21) preserves the invariant quadrants of the inviscid dynamics (just as in the case of uniform artificial viscosity) and provides the bound as long as for some . In particular, the Riemann invariants and satisfy and and the intersection of these half-planes provides the invariant quadrant (see Chueh, Conley, & Smoller ChCoSm1977 ()), and hence the desired bound as long as vacuum is avoided.

Finally, the fact that we have the lower-bound is an immediate consequence of the strong maximum principle.

### 2.10 The C-equation as a gradient flow

Notice that equilibrium solutions to the -equation are minimizers of the following functional (we drop the superscript ):

 EG(C)=∫(12C2x−G(ux)C+12ϵC2)dx.

In the absence of a forcing function , this reduces to

 E0(C)=12∫(C2x+1ϵC2)dx. (27)

The first term is commonly referred to as the Dirichlet energy and its minimizers are harmonic functions. The second term can be viewed as a penalization of the Dirichlet energy. In particular, because the energy functional is bounded by a constant independent of , the penalization term constrains to be . Thus, minimizers are trying to be harmonic while minimizing their support.

The -equation can be written as a classical gradient flow equation

 dCdt=−S(u)∇EG(C),

where the gradient is computed relative to the inner-product. Thus the heat operator in the -equation, , smooths the forcing in space-time, while the reaction term minimizes the support of the smoothed profile. This is very much related to the theories of Cahn-Hilliard and Ginzburg-Landau gradient flows, and we intend to examine this connection in subsequent papers.

## 3 Numerical Schemes

We describe two very different numerical algorithms in the context of our -method. First, we outline a classical continuous finite-element discretization, yielding FEM-C and FEM- (based on classical artificial viscosity). Second, we discuss a simple WENO-based scheme for compressible Euler that upwinds solely based on the sign of the velocity . To this scheme, we apply a slightly modified -method resulting in our WENO-C algorithm.

For the purpose of comparison, we also implement two additional numerical methods. The first is a second-order central-differencing scheme of Nessayhu-Tadmor (NT), a nice and simple method which serves as a base-line for our FEM-C comparisons. The second scheme is a very competitive WENO scheme that utilizes a Godunov-based upwinding based upon characteristic decompositions (WENO-G). This will serve as a benchmark for our WENO-C scheme.

### 3.1 Notation for discrete solutions

To compute approximations to (4), we subdivide space-time into a collection of spatial nodes and temporal nodes . We denote the computed approximate solution by

 uni≈u(xi,tn),

noting that for fixed and , is a 3-vector of solution components, i.e.,

 uni=⎡⎢ ⎢⎣ρnimniEni⎤⎥ ⎥⎦.

It is important to note that we use the notation for both pointwise approximations to , (acquired via FEM-C) and approximations to the cell-average values of (acquired via WENO-C).

A subscripted quantity denotes the vector itself and the individual components of the vector. We overload this notation so to not cause any confusion between functions defined over a continuum versus those defined only at a finite number of points.

In FEM-C and WENO-C, we discretize (9) (or some slight modification) with , and use the above notation for the computed solution. We also denote the approximation to by .

### 3.2 FEM-C and FEM-|ux|: A Second-Order Continuous-Galerkin Finite-Element Scheme

We choose a second-order continuous-Galerkin finite-element scheme to provide a discretization of (9), subsequently defining our FEM-C scheme.

We subdivide with (for even)-uniformly spaced nodes separated by a distance . In the FEM community, spatial discretization size is more commonly referred by element-width; to maintain consistency with the literature, we refer to the inter-nodal regions as cells. Since we use a continuous FEM, the degrees-of-freedom are defined at the cell-edges (as opposed to cell-centers)333When we compare our FEM-C scheme with other, cell-averaged schemes, we perform an averaging procedure based upon averages between nodes..

For use in our FEM implementation, it is useful to consider the variational form of (9). At the continuum level, satisfy

 ∫I⎡⎢⎣∂tuϵ⋅Φ−F(uϵ)⋅∂xΦ+βϵ2maxI |∂xuϵ|maxI CϵCϵ∂xuϵ⋅∂xΦ⎤⎥⎦dx=0, (28a) ∫I[∂tCϵϕ+S(uϵ)(ϵ∂xCϵ∂xϕ+1ϵCϵϕ)]dx=∫IS(uϵ)G(∂xuϵ)ϕ dx (28b)

for almost every , for all vector-valued test functions , and all scalar-valued test functions .

Using the finite-element spatial discretization based on piecewise second-order Lagrange polynomials, we construct operators and , corresponding to the non-time-differentiated terms in (28a) and (28b), respectively. Using these discrete operators, we write the semi-discrete form of (28a) and (28b) as

 ∂t[uiCi]+[AFEM(ui,Ci)BFEM(ui,Ci)]=0 (29)

where and represent the nodal values of an approximation to and for which (see Section 2.6). For a standard reference on the details of this procedure, see Larsson & Thomée Larsson2003 ().

The time-differentiation in (29) is approximated by a diagonally-implict second-order time-stepping procedure; first we predict to and solve the implicit set of equations for and follow by implicitly solving for using . Our fully discrete scheme is given by

 ~un+1i =uni+AFEM(uni,Cni), (30a) Cn+1i =Cni+tn+1−tn2[BFEM(~un+1i,Cn+1i)+BFEM(uni,Cni)], (30b) un+1i =uni+tn+1−tn2[AFEM(un+1i,Cn+1i)+AFEM(uni,Cni)]. (30c)

For smooth solutions, where artificial viscosity is not necessary, our FEM-C scheme is second-order accurate in both space and time when the error is measured in the -norm. Moreover, the addition the artificial viscosity obtained through the -method is formally a second-order perturbation (in ) and we have verified this accuracy when (again, for smooth ). For containing jump discontinuities, the given scheme is no longer second-order accurate on all of but preserves second-order accuracy in the smooth regions away from discontinuities.

For the classical artificial viscosity schemes (8), the C-equation is no longer used but we require a similar step to predict the velocity for use in the diffusion coefficient. This analogous scheme, is referred to as the FEM- scheme.

### 3.3 WENO-C: A Simple WENO scheme using the C-method

Our WENO-based scheme is motivated by Leonard’s finite volume schemes (Leonard (), pg. 65). Upwinding is performed solely based on the sign of the velocity at cell-edges, and the WENO reconstruction procedure is formally fifth-order.

We divide the interval into equally sized cells of width , identifying the degrees-of-freedom as cell-averages over cells centered at the . The cell edges are denoted using the fraction index, i.e.

 xi+1/2=xi+xi+12

Subsequently, we denote a cell-averaged quantity by and its values at the left and right endpoints by and , respectively.

Given a vector , corresponding to cell-average values, and vectors , corresponding to left and right cell-edge values, we define the th component of vector

 [WENO(wi,zi±1/2)]j=1Δx(~wj+1/2zj+1/2−~wj−1/2zj−1/2)

where the cell-edge values of are calculated using a fifth-order WENO reconstruction, upwinding based upon the sign of .

For the flux in the energy equation, we use

 [WENOE(Ei,ui±1/2)]j=1Δx(~Ej+1/2uj+1/2(1+pjEj)+(1+pj+1Ej+1)2−~Ej−1/2uj−1/2(1+pj−1Ej−1)+(1+pjEj)2). (31)

Using this simplified WENO-based reconstruction, we construct the operators and where

 ⎡⎢⎣AWENO⎛⎜⎝⎡⎢⎣ρimiEi⎤⎥⎦, Ci⎞⎟⎠⎤⎥⎦=⎡⎢ ⎢ ⎢⎣% WENO(ρi,ui±1/