On the scaling of entropy viscosity in high order methods

TL;DR

This paper examines how different scaling choices affect the entropy viscosity method's effectiveness in high order numerical schemes for shock problems, providing insights into optimal viscosity scaling strategies.

Contribution

It introduces a detailed analysis of entropy viscosity scaling and compares the performance of two scaling approaches through illustrative examples.

Findings

Scaling choice significantly impacts viscosity size and shock capturing quality

Different scalings lead to varied numerical stability and accuracy

Guidelines for selecting appropriate entropy viscosity scaling are discussed

Abstract

In this work, we outline the entropy viscosity method and discuss how the choice of scaling influences the size of viscosity for a simple shock problem. We present examples to illustrate the performance of the entropy viscosity method under two distinct scalings.

Tables

Table 1. Table 1: Size of μ E subscript 𝜇 𝐸 \mu_{E} and μ max subscript 𝜇 max \mu_{\rm max} for different size of shock speed ( v s subscript 𝑣 𝑠 v_{s} ) with respect to the size of the jump ( Δ ​ u Δ 𝑢 \Delta u ) in the entropy viscosity method.

Case	${\nu }_{EV}$	${\nu }_{\max }$
\svhline $\left|v_S\right|\ll \left|\Delta u\right|$	${\alpha }_{EV}{h}^{\beta -1}{\left|\Delta u\right|}^3$	${\alpha }_{\max }h\left|v_s\right|$
$\left|v_s\right|\approx \left|\Delta u\right|$	${\alpha }_{EV}{h}^{\beta -1}{\left|\Delta u\right|}^3$	$2{\alpha }_{\max }h\left|v_s\right|$
$\left|v_s\right|\gg \left|\Delta u\right|$	${\alpha }_{EV}{h}^{\beta -1}{\left|\Delta u\right|}^3$	$0.5{\alpha }_{\max }h\left|\Delta u\right|$

Equations

$$u_t+f(u{)}_x=0,$$

$$r_{EV}=E_t+F_x\equiv {\left(\frac{u^2}{2}\right)}_t+{\left(\frac{u^3}{3}\right)}_x\le 0.$$

$$\nu =\min ({\nu }_{\max },{\nu }_{EV}),$$

$$u_t+f(u{)}_x=(\nu u_x{)}_x.$$

$${\nu }_{EV}(x)={\alpha }_{EV}{h}^{\beta }|r_{EV}(x)|,$$

$$\frac{d{u}_i}{dt}\approx \frac{f_{i+\frac{1}{2}}-f_{i-\frac{1}{2}}}{h},$$

$$f_{i+\frac{1}{2}}(u_i,u_{i+1})=\frac{1}{2}{\left(\frac{u_i+u_{i+1}}{2}\right)}^2-\left(\frac{{\nu }_i+{\nu }_{i+1}}{2}\right)\frac{u_{i+1}-u_i}{h}.$$

$$\frac{d{F}_i}{dx}=\frac{F_{i+\frac{1}{2}}-F_{i-\frac{1}{2}}}{h},\text{ where, }{F}_{i+\frac{1}{2}}=\frac{1}{3}{\left(\frac{u_i+u_{i+1}}{2}\right)}^3.$$

$${\nu }_{\max }={\alpha }_{\max }h\max |u|,$$

$$\nu=\left\{\begin{array}[]{ll}0&s<s_{0}-\kappa,\\ \varepsilon_{0}h&s>s_{0}+\kappa,\\ \frac{\varepsilon_{0}h}{2}\left(1+\sin\left(\frac{\pi(s-s_{0})}{2\kappa}\right)\right)&\text{otherwise}.\end{array}\right.$$

$$u(x,t)=u_l+\Delta uH\left(x-v_{s}t\right).$$

$$u_t+{\left(\frac{u^2}{2}\right)}_x$$

$$=-\Delta u\left(\frac{2u_l+\Delta u}{2}\right)H^{\prime }+\Delta u\left(\frac{2u_l+\Delta u}{2}\right)H^{\prime }$$

$$=0,$$

$$r_{EV}=\frac{(\Delta u{)}^3}{12}{H}^{\prime }(x-x_s)=\frac{(\Delta u{)}^3}{12}\delta (x-x_s).$$

$${\int }_{-\infty }^{\infty }\delta (x)dx=1.$$

$$\sum _{j=0}^{n-1}{h}_j{\delta }_j=1,$$

$${\nu }_{EV}={\alpha }_{EV}{h}^{\beta -1}|(\Delta u{)}^3|\text{ and }{\nu }_{\max }={\alpha }_{\max }h\max (|u_l|,|u_r|).$$

$$u(x,t)=\{\begin{array}{ll}-0.5+v_s,& x\in [-1,v_{s}t),\\ 0.5+v_s,& x\in [v_{s}t,1],\end{array}$$

$$u(x,0)=-\sin (\pi x)+0.5,$$

$$u(x,0)=\{\begin{array}{ll}0& x\in [-1,-0.5),\\ -0.1\sin (2\pi x)& x\in [-0.5,0.5),\\ 0& x\in [0.5,4.5),\\ -0.5& x\in [4.5,5],\end{array}$$

Full text

1 Introduction

Hyperbolic partial differential equations (PDE) are used to model various fluid flow problems. In the special case of 1-dimensional linear constant coefficient scalar hyperbolic problems, the solutions to these PDE are simply a translation of the initial data. However, for nonlinear problems the solution may deform, and as a result, shock waves can form even if the initial data is smooth leveque2002finite .

In computational fluid dynamics, it is desirable that numerical methods capture shock waves and maintain a high accuracy for smooth waves. Low order methods have sufficient numerical dissipation to regularize shock waves but obtaining accurate solutions in smooth regions can become expensive. On the other hand, high order methods are capable of achieving high accuracy at a reasonable cost. Their low numerical dissipation enables such accuracy, but on the downside, it limits their ability to regularize shock waves.

Various techniques have been implemented to capture shocks while maintaining high accuracy, at least away from shocks. There are two major classes of shock capturing techniques: shock detection techniques, where we find slope limiters leveque2002finite , Essentially Non-Oscillatory (ENO) and Weighted ENO (WENO) shu1998essentially , and artificial viscosity techniques, where we find filtering yee2007development ; persson2006sub , the PDE-based viscosity method johnson1990convergence , the entropy viscosity method Guermond2011conservation , among others.

In this work, we focus on the entropy viscosity method. In essence, the entropy viscosity method provides shock capturing without compromising the high accuracy away from the shock. An important advantage of this method is that it generalizes very easily to higher dimensions and unstructured grids.

As a model problem, we consider Burgers’ equation

[TABLE]

where $f=\frac{u^{2}}{2}$. Physically correct solutions to (1) can be singled out by requiring that they satisfy an entropy inequality such as

[TABLE]

The entropy residual, $r_{EV}$, is zero wherever $u$ is smooth. If the solution $u$ contains a shock, then the entropy residual takes the form of a negative Dirac distribution centered at the location of the shock, $x_{s}$, i.e. $r_{EV}=-C\,\delta(x-x_{s})$. The property that the entropy residual is unbounded at a shock was first used by Guermond and Pasquetti in Guermond2011conservation , as a way to selectively introduce viscosity. The artificial viscosity, $\nu$, proposed in Guermond2011conservation , defined as the minima of two viscosities

[TABLE]

becomes the coefficient of the viscous term in the viscous Burgers’ equation,

[TABLE]

Here, $\nu_{\rm max}$ is the Lax-Friedrich viscosity whose size depends on discretization and the largest eigenvalue, $\lambda_{\rm LF}$, of the flux Jacobian, $\frac{Df(u)}{Du}$. The second viscosity $\nu_{EV}$ is proportional to the magnitude of the entropy residual (in fact, a discretization of the entropy residual) and will thus be zero (or small after discretization) away from discontinuities. In theory, the entropy residual becomes unbounded at a shock, numerically however, the entropy residual $r_{EV}$ remains bounded with the size of the residual depending on the discretization size. As we will see below, this subtle difference has consequences for how to choose the scaling of the viscosity terms in the entropy viscosity method.

On a grid with step size $h$, the second viscosity $\nu_{EV}$ can be expressed as

[TABLE]

with a parameter $\alpha_{EV}$ that requires tuning. In recent papers on entropy viscosity method, see e.g. guermond2011entropyhighorder ; guermond2010entropy ; guermond2011entropy ; guermond2011suitable ; zingan2013implementation , the parameter $\beta$ is chosen to be 2, but the original paper guermond2008entropy uses $\beta=1$. It is unclear to us why the later works prefer $\beta=2$. Here, we will present analysis and computational results that suggest the original scaling $\beta=1$ is a more natural choice. We note that the entropy residual is typically scaled by $\|E-\overline{E}\|_{\infty}$, with the over-bar indicating a spatial average, but as this quantity is roughly constant in the problems presented here, we omit it for brevity and reduced complexity.

The rest of the paper is organized as follows. In Section 2, we describe different discretizations of (4) that we consider here, in Section 3, we present an analysis of how the entropy viscosity $\nu$ depends on the two viscosities, $\nu_{EV}$ and $\nu_{\rm max}$, under different scaling for a model problem. In Section 4, we then conduct experiments with the entropy viscosity method where $\beta$ takes on values 1 or 2 and compare the results.

2 Numerical methods

We will consider the discretization of (4) by our conservative Hermite method KorAppJSC , a standard discontinuous Galerkin (dG) method Hesthaven:2008fk and a simple finite volume type discretization leveque2002finite . For all the discretizations we let the domain $x_{L}\leq x\leq x_{R}$ be discretized by the regular grid $x_{i}=x_{L}+i\,h,\,i=0,\ldots,n,$ $h=(x_{R}-x_{L})/n$.

The degrees of freedom for the finite volume method is cell averages centered at the grid points. For the Hermite method, the degrees of freedom are the coefficients of node centered Taylor polynomials of degree $m$ and for the dG method, they are the ($m+1$) coefficients of element-wise (we take an element to be $\Omega_{i}=[x_{i-1},x_{i}]$) expansions in Legendre polynomials. For smooth solutions the spatial accuracy of the Hermite method is $2m+1$ and $m+1$ for the dG method.

All three methods use the classic fourth order Runge-Kutta method to evolve the semi-discretizations in a method-of-lines fashion.

In the Hermite method, we evaluate the fluxes and their derivatives at the nodes (element edges) for the four stages in the RK method. Precisely, for the first stage we compute the slope $f^{h}_{1}=\frac{1}{2}\mathcal{T}[(u^{h}_{1})^{2}]-\frac{\nu}{h}\frac{du^{h}_{1}}{dx}$ for the Taylor polynomial $u^{h}_{1}=u^{h}$ approximating the solution at the first stage. Here $\mathcal{T}[(u^{h}_{1})^{2}]$ is the truncated polynomial multiplication of $u^{h}_{1}$ with itself and $\frac{du^{h}_{1}}{dx}$ is the derivative of the polynomial. At the next stage, the solution is $u^{h}_{2}=u^{h}+\frac{(\Delta t/2)}{2}\frac{df^{h}_{1}}{dx}$, the slope is $f^{h}_{2}=\frac{1}{2}\mathcal{T}[(u^{h}_{2})^{2}]-\frac{\nu}{h}\frac{du^{h}_{2}}{dx}$ and so on. Once the stage slopes $f_{s}^{h},\,s=1,\ldots,4$ and their spatial derivatives are known, we perform a Hermite interpolation to the element centers of the solution and the spatial derivatives of the stage slopes. These are then used to evolve the element centered Hermite interpolant of $u^{h}$ to $t=t_{n}+\Delta t/2$. As the Hermite interpolant is of higher degree than the original Taylor polynomial, we conclude a half-step by truncating it to the appropriate degree. To advance the solution a full time step, the half-step process is repeated starting from the element centers.

To handle the artificial viscosity in the dG method, we use the approach of Bassi and Rebay bassi1997high with a Lax-Friedrichs flux for the advective term and alternating fluxes for the viscous term. The nonlinear terms are constructed explicitly and de-aliased by over-integration kirby2003aliasing .

For the finite volume method, we let $u_{i}\approx u(x_{i})$ be a grid function approximating the solution and $f_{i+\frac{1}{2}}=f_{i+\frac{1}{2}}(u_{i},u_{i+1})$ be an approximation to the flux at $x_{i+\frac{1}{2}}$. To compute the time derivatives, we use the spatial approximation

[TABLE]

where

[TABLE]

When $\nu_{i}=0$, the above discretization is linearly stable (when paired with a suitable time-stepping method) but is not non-linearly stable, and we thus add artificial viscosity to stabilize it.

For all three discretizations, we approximate the time derivative of the entropy function, $E_{t}$, by a backward difference. This approach is explicit as we use the current solution to compute $E$ at the current time before evolving the solution in time. The residual (and hence the viscosity) is kept on each element / grid-point over each step.

To approximate the entropy flux derivative $F_{x}$ using the Hermite method, we compute the derivative of the truncated polynomial multiplication $\mathcal{T}[u^{h}\mathcal{T}[(u^{h})^{2}]]$ at the node. For the dG method, we evaluate the flux $F$ on a Legendre-Gauss-Lobatto (LGL) grid and differentiate it to get an approximation for $F_{x}$. The residual on an element is taken to be the maximum of the absolute value of the residual on the LGL grid. In the finite volume method $F_{x}$ is approximated by

[TABLE]

We note that more sophisticated discretizations of the entropy residual could be considered. In particular, a higher order approximation to $r_{EV}$ would result in a higher rate of convergence for smooth solutions, but as we are mainly concerned with the scaling $\beta$, we did not pursue such discretizations here. In fact, in our experience, the results concerning the choice of scaling are not affected by the order of the accuracy of the approximation to $r_{EV}$. This will be discussed in Section 3.

We also define $\nu_{\rm max}$ to be the classical Lax-Friedrich viscosity, which for Burgers’ equation takes the form

[TABLE]

where the maximum is taken globally.

Finally, for the purpose of comparison we also present some results computed using the sub-cell resolution smoothness sensor of Persson and Peraire, persson2006sub . The smoothness sensor compares the $L_{2}$ energy content of the highest (Fourier or expansion) mode with the total $L_{2}$ energy on an element and then maps its ratio (which is an indicator of the smoothness) into the size of the artificial viscosity. Precisely, if the approximate dG solution on an element is $u^{h}=\sum_{k=0}^{m}\hat{u}_{k}P_{k}$, with $P_{k}$ being an orthogonal basis, we compute the smoothness as $s=\log_{10}(\|\hat{u}_{m}P_{m}\|^{2}/\|u^{h}\|^{2})$ and the viscosity as

[TABLE]

When applied to the Hermite method, we first project the Taylor polynomials centered at two adjacent grid-points into an orthogonal Legendre expansion on the element defined by the grid-points and then proceed as above.

3 Impact of the $h$-scaling on the selection mechanism

To study how the selection mechanism depends on the shock speed and the size of the jump, consider a solution of the Burgers’ equation consisting of a Heaviside function $H$ with left state $u_{l}$ and right state $u_{r}$, given by

[TABLE]

This corresponds to a shock of size $|\Delta u|=|u_{r}-u_{l}|$ moving with speed $v_{s}=0.5(u_{l}+u_{r})$. Solutions of the form (9) always has a negative $\Delta u$ value since Lax entropy condition for Burgers’ equation dictates $u_{l}=f^{\prime}(u_{l})>v_{s}>f^{\prime}(u_{r})=u_{r}$.

For simplicity, we use the short hand notation $H$ for $H\left(x-v_{s}t\right)$. A direct computation

[TABLE]

shows that (9) is a solution of (1). Further, it can be shown that the entropy residual (2) for (9) is

[TABLE]

That is, the size of the entropy residual grows with the cube of $\Delta u$.

Now, by the properties that define the Dirac delta function $\delta$, we have

[TABLE]

Thus, a consistent discretization of the Dirac delta function $\delta_{0},...,\delta_{n}$ on a grid $x_{0},...,x_{n}$ must obey the condition

[TABLE]

where $h_{j}=x_{j+1}-x_{j}$. For any approximation with a finite width stencil, we must have $\delta_{j}\sim h^{-1}_{j}$ and we thus expect $r_{EV}$ to behave like $(\Delta u)^{3}/h$ on a uniform grid. We therefore proceed with the analysis using the discrete approximation $r_{EV}=(\Delta u)^{3}/h$. Using this expression for $r_{EV}$, we estimate the viscosity $\nu$ by the minimum of

[TABLE]

The comparison between the size of $\nu_{EV}$ and $\nu_{\rm max}$ in various scenarios is reported in Table 3. If $\beta=2$, then the two viscosities $\nu_{EV}$ and $\nu_{\rm max}$ scale as $h$. For a problem with multiple shocks, the homogeneity in $h$-scaling introduces an additional difficulty in determining ${\alpha}_{EV}$. Should it be chosen based on the largest or smallest shock? What if new shocks appear during the course of the computation? To avoid answering these questions, we instead consider $\beta=1$. Now $\nu_{EV}=\mathcal{O}(1)$ while $\nu_{\rm max}=\mathcal{O}(h)$, and the particular choice of ${\alpha}_{EV}$ is thus irrelevant since as $h\to 0$, the selection mechanism will eventually select $\nu_{\rm max}$ at the shocks. We will provide an example to illustrate the two-shock dilemma in Section 4.3.