Shock Capturing by Bernstein Polynomials for Scalar Conservation Laws

TL;DR

This paper introduces a novel shock capturing method using Bernstein polynomials to stabilize high-order spectral element solutions of hyperbolic conservation laws, effectively reducing oscillations near discontinuities without compromising computational efficiency.

Contribution

The authors develop a Bernstein-based shock capturing procedure that is total variation diminishing, preserves monotonicity, and can enforce bounds, improving stability of spectral methods for conservation laws.

Findings

The method stabilizes spectral approximations near shocks.

It preserves monotone and bounded solutions.

It does not reduce the time step or CFL condition.

Abstract

A main disadvantage of many high-order methods for hyperbolic conservation laws lies in the famous Gibbs-Wilbraham phenomenon, once discontinuities appear in the solution. Due to the Gibbs-Wilbraham phenomenon, the numerical approximation will be polluted by spurious oscillations, which produce unphysical numerical solutions and might finally blow up the computation. In this work, we propose a new shock capturing procedure to stabilise high-order spectral element approximations. The procedure consists of going over from the original (polluted) approximation to a convex combination of the original approximation and its Bernstein reconstruction, yielding a stabilised approximation. The coefficient in the convex combination, and therefore the procedure, is steered by a discontinuity sensor and is only activated in troubled elements. Building up on classical Bernstein operators, we are thus able to prove that the resulting Bernstein procedure is total variation diminishing and preserves monotone (shock) profiles. Further, the procedure can be modified to not just preserve but also to enforce certain bounds for the solution, such as positivity. In contrast to other shock capturing methods, e.g. artificial viscosity methods, the new procedure does not reduce the time step or CFL condition and can be easily and efficiently implemented into any existing code. Numerical tests demonstrate that the proposed shock-capturing procedure is able to stabilise and enhance spectral element approximations in the presence of shocks.

Equations

$${\partial }_{t}u+{\partial }_{x}f(u)=0,x\in \Omega \subset \mathbb{R},$$

$${\partial }_{t}u+a(u){\partial }_{x}u=0$$

$${\partial }_{t}U(u)+{\partial }_{x}F(u)\le 0$$

$$b_{n,N}(x)=\left(\genfrac{}{}{0ex}{}{N}{n}\right)\frac{(x-a{)}^n(b-x{)}^{N-n}}{(b-a{)}^n},n=0,\dots ,N,$$

$$B_N(x)=\sum _{n=0}^N{\beta }_n{b}_{n,N}(x),$$

$$B_N[u](x)=\sum _{n=0}^{N}u\left(a+\frac{n}{N}(b-a)\right)b_{n,N}(x).$$

$${\beta }_n=u\left(a+\frac{n}{N}(b-a)\right).$$

$$\sum _{n=0}^N{b}_{n,N}(x)=1$$

$$m\le {\beta }_n\le M\forall n=0,\dots ,N⟹m\le B_N(x)\le M$$

$$m\sum _{n=0}^N{b}_{n,N}(x)\le \sum _{n=0}^N{\beta }_n{b}_{n,N}(x)\le M\sum _{n=0}^N{b}_{n,N}(x)$$

$$B_N[u](0)=u(0)\text{and}{B}_N[u](1)=u(1)$$

$$B_N^{(k)}[u](x)=N(N-1)\dots (N-k+1)\sum _{n=0}^{N-k}{\Delta }^{k}u\left(\frac{n}{N}\right)\left(\genfrac{}{}{0ex}{}{N-k}{n}\right)x^n(1-x{)}^{N-n-k}$$

$${\Delta }^{k}u\left(\frac{n}{N}\right)=\Delta \left({\Delta }^{k-1}u\left(\frac{n}{N}\right)\right)=u\left(\frac{n+k}{N}\right)-\left(\genfrac{}{}{0ex}{}{k}{1}\right)u\left(\frac{n+k-1}{N}\right)+\cdots +(-1{)}^{k}u\left(\frac{n}{N}\right)$$

$$B_N^{\prime }[u](x)=N\sum _{n=0}^{N-1}\left[u\left(\frac{n+1}{N}\right)-u\left(\frac{n}{N}\right)\right]\left(\genfrac{}{}{0ex}{}{N-1}{n}\right)x^n(1-x{)}^{N-1-n}$$

$$\underset{N\to \infty }{\lim }N\left(u(x)-B_N[u](x)\right)=-\frac{x(1-x)}{2}{u}^{\prime \prime }(x)\text{for}x\in [0,1];$$

$$\left|u(x)-B_N[u](x)\right|\le C{N}^{-1}{u}^{\prime \prime }(x)$$

$$u(x_k^{+}):=\underset{x\to x_k^{+}}{\lim }u(x)\text{and}u(x_k^{-}):=\underset{x\to x_k^{-}}{\lim }u(x)$$

$${\left({\int }_0^1{\left|u(x)-B_N[u](x)\right|}^p\mathrm{d}x\right)}^{\frac{1}{p}}\le C{N}^{-1}$$

$${\int }_0^1{\left|u(x)-B_N[u](x)\right|}^p\mathrm{d}x=\sum _{k=0}^K{\int }_{x_k}^{x_{k+1}}{\left|u(x)-B_N[u](x)\right|}^p\mathrm{d}x$$

$$\left|u(x)-B_N[u](x)\right|\le C{N}^{-1}{u}^{\prime \prime }(x)$$

$$\left|u(x)-B_N[u](x)\right|\le C{N}^{-1}$$

$${\int }_{x_k}^{x_{k+1}}{\left|u(x)-B_N[u](x)\right|}^p\mathrm{d}x$$

$$={\int }_{x_k}^{x_k+\epsilon }{\left|u(x)-B_N[u](x)\right|}^p\mathrm{d}x$$

$$+{\int }_{x_k+\epsilon }^{x_{k+1}-\epsilon }{\left|u(x)-B_N[u](x)\right|}^p\mathrm{d}x+{\int }_{x_{k+1}-\epsilon }^{x_{k+1}}{\left|u(x)-B_N[u](x)\right|}^p\mathrm{d}x$$

$$={\int }_{x_k}^{x_k+\epsilon }{\left|u(x)-B_N[u](x)\right|}^p\mathrm{d}x+(x_{k+1}-x_k-2\epsilon )C^p{N}^{-p}+{\int }_{x_{k+1}-\epsilon }^{x_{k+1}}{\left|u(x)-B_N[u](x)\right|}^p\mathrm{d}x.$$

$$\left|u(x)-B_N[u](x)\right|\le \left|u(x)\right|+\left|B_N[u](x)\right|\le 2{\Vert u\Vert }_{\infty }$$

$${\int }_{x_k}^{x_{k+1}}{\left|u(x)-B_N[u](x)\right|}^p\mathrm{d}x\le (x_{k+1}-x_k-2\epsilon )C^p{N}^{-p}+4\epsilon {\Vert u\Vert }_{\infty }$$

$${\int }_{x_k}^{x_{k+1}}{\left|u(x)-B_N[u](x)\right|}^p\mathrm{d}x\le (x_{k+1}-x_k)C^p{N}^{-p}.$$

$${\int }_{x_k}^{x_{k+1}}{\left|u(x)-B_N[u](x)\right|}^p\mathrm{d}x$$

$$⟹$$

$$B_N[u](x)=\sum _{n=0}^{N}u\left(-1+2\frac{n}{N}\right)b_{n,N}(x).$$

$${\underline{u}}^{(B)}=\underline{\underline{T}}\underline{b},$$

Full text

1 Introduction

In this work, we introduce a shock capturing procedure for spectral element (SE) approximations of scalar hyperbolic conservation laws

[TABLE]

with appropriate initial and boundary conditions. Here, $u(t,x)$ is the unknown function and $f(u)$ is a $C^{1}$ flux function. We assume that the initial data $u_{0}$ is a function of bounded total variation in $\Omega$. Equation (1) might as well be written in its quasi-linear form

[TABLE]

with $a(u)=f^{\prime}(u)$. We further recall that solutions of (1) may develop spontaneous jump discontinuities (shock waves and contact discontinuities) even when the initial data are smooth [35]. This important discovery has first been made by Riemann [50]. Hence, the more general class of weak solutions is admitted, where (1) is satisfied in the sense of distribution theory. Since there are many possible weak solutions, however, equation (1) is augmented with an additional entropy condition, requiring

[TABLE]

to hold. $U$ is an entropy function and $F$ is a corresponding entropy flux satisfying $U^{\prime}f^{\prime}=F^{\prime}$. A strict inequality in (3) reflects the existence of physically reasonable shock waves.

Next, we consider spectral/$hp$ element approximations of (1). Also see [29] and references therein. SE approximations have been introduced in 1984 by Patera [42] and can be interpreted as a formulation of finite elements (FE) that uses piecewise polynomials of high degrees. SE methods thus combine the advantages of (Galerkin) spectral methods with those of FEs by a simple application of the spectral method to subdomains (elements) $\Omega_{i}$ of $\Omega$. The elements $\Omega_{i}$ are mapped to a reference element $\Omega_{\mathrm{ref}}$ and all computations are carried out there. In one space dimension, this reference element is typically given by $\Omega_{\mathrm{ref}}=[-1,1]$. Thus, the SE approximation is based on a polynomial approximation of degree at most $N$, expressed w. r. t. some nodal or modal basis of $\mathbb{P}_{N}(\Omega_{\mathrm{ref}})$. This polynomial approximation, in particular, is used to extrapolate the solution to the element boundaries (if these values are not already given as coefficients in a nodal basis). Subsequently, common numerical fluxes $f^{\mathrm{num}}$ are computed at each element boundary [55]. Numerical fluxes allow information to cross element boundaries and are often essential for the stability of the method. They are incorporated into the approximation by some correction procedure at the boundaries. In discontinuous Galerkin (DG) methods, this is done by applying integration by parts to a weak form of (1), see [25]. In flux reconstruction (FR) methods, which are based on the strong form of (1), this is done by applying correction functions, see [28]. Finally, $\partial_{x}f(u)$ is approximated using exact differentiation for the polynomial approximation.

Convergence of the resulting SE method is either obtained by increasing the degree of the polynomials ($p$-refinement) or by increasing the number of elements ($h$-refinement). SE methods are known to provide highly accurate approximations in smooth regions. But also from a computational point of view, the main difficulty of solving hyperbolic conservation laws is that discontinuities may arise: Resulting from the famous Gibbs–Wilbraham phenomenon [26], the polynomial approximation is polluted by spurious oscillations, which might produce unphysical numerical solutions and finally blow up the computation. To overcome these problems, many researchers have extended shock capturing procedures from finite difference (FD) as well as finite volume (FV) methods to (high-order) SE schemes. The idea behind many of these procedures is to add dissipation to the numerical solution. This idea dates back to the pioneering work [57] of von Neumann and Richtmyer during the Manhattan project in the 1940’s at Los Almos National Laboratory, where they constructed stable FD schemes for the equations of hydrodynamics by including artificial viscosity (AV) terms. Since then, AV methods have attracted the interest of many researchers and were investigated also for SE approximations in a large number of works [43, 31, 19, 48, 18]. Despite providing a robust and accurate way to capture (shock) discontinuities, AV terms are not trivial to include in SE approximations. Typically, they are nonlinear and consist of higher (second and fourth order) derivatives. Another drawback arises from the fact that AV terms can introduce additional harsh time step restrictions, when not constructed with care, and thus decrease the efficiency of the numerical method [25, 18]. Finally, we mention those methods based on order reduction [5, 10], mesh adaptation [14], weighted essentially non-oscillatory (WENO) concepts [52, 53], and $\ell^{1}$ regularisation applied to high order approximations of the jump function [17]. Yet, a number of issues still remains unresolved. Often, the extension of these methods to multiple dimensions is an open question or they are regarded as too computational expensive.

Here, we propose an approach that combines the good properties of SE approximations in smooth regions with the total variation diminishing and shape preserving properties of Bernstein approximations for resolving shocks without spurious oscillations. In this new strategy, the approximation of $u$ in each element may vary from usual (high-order) interpolation polynomials to a Bernstein reconstruction of the solution. Further, by employing a discontinuity sensor, here based on comparing polynomial annihilation operators [3] of increasing orders as proposed in [17], the order of the approximations is reduced to one only in elements where the solution is not smooth. For instance mesh adaptation is hence not mandatory and (shock) discontinuities can be captured without modifying the number of degrees of freedom, the mesh topology, or even the method. Moreover, a slight modification of the proposed Bernstein reconstruction also allows the preservation of bounds, such as positivity of pressure and water height.

By now, Bernstein polynomials as polynomial bases have been successfully applied to solvers for partial differential equations (PDEs) in a number of works. Especially in [30, 1], efficient FE operators have been constructed by using Bernstein polynomials as shape functions. Also in [6], Bernstein polynomials have been proposed as basis functions in a third-order quasi-nodal DG method for solving flooding and drying problems with shallow water equations. Finally, [38, 2] have investigated the potential of imposing discrete maximum principles for polynomials expressed in a basis of Bernstein polynomials. Yet, while Bernstein polynomials as basis functions have been studied in a few works by now, the associated approximation procedure resulting from the Bernstein operator is still of very limited use. In contrast to the above mentioned works, we do not only utilise Bernstein polynomials as basis functions, but we further propose to use the associated Bernstein operatore as a building stone for new shock capturing procedures.

The rest of this work is organised as follows. In §2, we revise bases of Bernstein polynomials and their associated Bernstein approximation operator. This operator will later be used to obtain ’smoother’ reconstructions of polynomial approximations near discontinuities and provides certain structure-preserving and approximation properties. Building up on these properties, §3 introduces the novel Bernstein procedure, which replaces polluted interpolation polynomial by convex combinations of the original (polluted) approximation and its Bernstein reconstruction. This process is steered by a polynomial annihilation sensor, which is discussed in §3.4. §4 investigates some analytical properties of the proposed Bernstein procedure, such as its effect on entropy, total variation and monotone (shock) profiles of the numerical solution. We stress that the Bernstein reconstruction is proven to be total variation diminishing (nonincreasing). Finally, §5 provides numerical demonstrations for a series of different scalar test problems. We close this work with concluding thoughts in §6.

2 Bernstein Polynomials and the Bernstein Operator

In this section, we introduce Bernstein polynomials as well as some of their more important properties. On an interval $[a,b]$, the $N+1$ Bernstein basis polynomials of degree $N$ are defined as

[TABLE]

and form a basis of $\mathbb{P}_{N}$. Thus, every polynomial of degree at most $N$ can be written as a linear combination of Bernstein basis polynomials,

[TABLE]

called Bernstein polynomial or polynomial in Bernstein form of degree $N$. The coefficients are referred to as Bernstein coefficients or Bézier coefficients. We further define the linear Bernstein operator of degree $N$ for a function $u\mathrel{\mathop{\ordinarycolon}}[a,b]\to\mathbb{R}$ by

[TABLE]

$B_{N}[\cdot]$ maps a function $u$ to a Bernstein polynomial of degree $N$ with Bernstein coefficients

[TABLE]

Without loss of generality, we can restrict ourselves to the intervals $[0,1]$ and $[-1,1]$. For sake of simplicity, the interval $[0,1]$ subsequently will be used for theoretical investigations. The interval $[-1,1]$, on the other hand, is typically used as a reference element in SE approximations. Thus, the proposed Bernstein procedure will be explained for this case.

2.1 Structure-Preserving Properties

Let us consider Bernstein polynomials on $[0,1]$. We start by noting that the Bernstein basis polynomials form a partition of unity, i. e.

[TABLE]

for all $N\in\mathbb{N}$. This is a direct consequence of the binomial theorem. Thus, we can immediately note

Lemma 1**.**

Let $B_{N}(x)=\sum_{n=0}^{N}\beta_{n}b_{n,N}(x)$ be a Bernstein polynomial of degree $N$ with Bernstein coefficients $\beta_{1},\dots,\beta_{N}$. Then

[TABLE]

holds for the Bernstein polynomial.

Proof.

Let $m\leq\beta_{n}\leq M$ for all $n=0,\dots,N$. We therefore have

[TABLE]

and the assertion follows from (8). ∎

In particular, Lemma 1 ensures that the Bernstein operator (6) preserves the bounds of the underlying function $u$. In fact, Lemma 1 not only ensures preservation of bounds by the Bernstein procedure, but also allows us to enforce such bounds. Moreover, the Bernstein operator preserves the boundary values of $u$, i. e.

[TABLE]

hold. This makes the later proposed Bernstein procedure a reasonable shock-capturing method not just in discontinuous FE approximations, such as DG methods, but also in continuous FE approximations, where the numerical solution is required to be continuous across element interfaces. Further, we revise the formula

[TABLE]

with forward difference operator

[TABLE]

for derivatives of (6); see [39, Chapter 1.4]. In particular, we have

[TABLE]

for the first derivative of $B_{N}[u]$. This formula will be important later in order to show that the Bernstein procedure is able to preserve monotone (shock) profiles.

2.2 Approximation Properties

Bernstein polynomials were first introduced by Bernstein [8] in a constructive proof of the famous Weierstrass theorem, which states that every continuous function on a compact interval can be approximated arbitrarily accurate by polynomials [58]. Hence, the sequence of Bernstein polynomials $\left(B_{N}[u]\right)_{N\in\mathbb{N}}$ converges uniformly to the continuous function $u$. Assuming that $u$ is bounded in $[0,1]$ and that the second derivative $u^{\prime\prime}$ exists at the point $x\in[0,1]$, we have

[TABLE]

see [39, chapter 1.6.1]. In particular, at points $x$ where the second derivative exists, the error of the Bernstein polynomial $B_{N}[u]$ therefore is of first order, i. e.

[TABLE]

for a $C>0$. However, we should remember that solutions of scalar hyperbolic conservation laws – which we intend to approximate – might contain discontinuities. The structure of these solutions has been determined by Oleinik [41, 40], Lax [34], Dafermos [13], Schaeffer [51], Tadmor and Tassa [54], and many more. Most notably for our purpose, Tadmor and Tassa [54] showed for scalar convex conservation laws that if the initial speed has a finite number of decreasing inflection points then it bounds the number of future shock discontinuities. Thus, in most cases the solution consists of a finite number of smooth pieces, each of which is as smooth as the initial data. Note that it is this type of regularity which is often assumed – sometimes implicitly – in the numerical treatment of hyperbolic conservation laws. Hence, it appears to be more reasonable to investigate convergence of the sequence of Bernstein polynomials $\left(B_{N}[u]\right)_{N\in\mathbb{N}}$ for only piecewise smooth functions $u$. Here, we call a function $u$ piecewise $C^{k}$ on $[0,1]$ if there is a finite set of points ${0<x_{1}<\dots<x_{K}<1}$ such that $\left.u\right|_{[0,1]\setminus\{x_{1},\dots,x_{K}\}}\in C^{k}$ and if the one-sided limits

[TABLE]

exist for all $k=1,\dots,K$. In this case $u$ is still integrable and the first order convergence of the Bernstein polynomials carries over in an $L^{p}$-sense.

Theorem 2**.**

Let $u$ be piecewise $C^{2}$ and let $\left(B_{N}[u]\right)_{N\in\mathbb{N}}$ be the sequence of corresponding Bernstein polynomials. Then there is a constant $C>0$ such that

[TABLE]

holds for all $N\in\mathbb{N}$ and $1\leq p<\infty$.

Proof.

Let $0<x_{1}<\dots<x_{k}<1$ be the points where $u$ is not $C^{2}$. When further denoting $x_{0}=0$ and $x_{K+1}=1$, we have

[TABLE]

and, since $u$ is piecewise $C^{2}$, there is a generic constant $C>0$ such that

[TABLE]

holds for $x_{k}<x<x_{k+1}$. It should be noted however that $u^{\prime\prime}$ might be unbounded on the open interval $(x_{k},x_{k+1})$. Therefore, let us choose $\varepsilon>0$ such that $\varepsilon<\min_{k=0,\dots,K}|x_{k}-x_{k+1}|$ and let us consider the closed interval $[x_{k}+\varepsilon,x_{k+1}-\varepsilon]$. On these intervals $u^{\prime\prime}$ is bounded and (20) becomes

[TABLE]

for $x_{k}+\varepsilon\leq x\leq x_{k+1}-\varepsilon$. Hence, we have

[TABLE]

For the two remaining integrals, we remember that $u$ and therefore all $B_{N}[u]$ are bounded. This yields

[TABLE]

for all $x\in[0,1]$ and $N\in\mathbb{N}$. As a consequence, (22) reduces to

[TABLE]

and letting $\varepsilon\to 0$ results in

[TABLE]

Finally, we have

[TABLE]

which yields the assertion. ∎

One might reproach that this is an unacceptable order of convergence. We reply to this argument with a few selected counter-arguments.

Remark 3**.**

• •

For numerical solutions of hyperbolic conservation laws, it is almost universally accepted that near shocks, the solution can be first order accurate at most [43]. Thus, accuracy of the numerical solution won’t decrease noticeably by reconstructing it as a Bernstein polynomial.

• •

In fact, high-order methods often not even provide accuracy of first order but constant (or even decreasing) accuracy. This is due to the Gibbs–Wilbraham phenomenon [26, 49, 20] for (polynomial) higher-order approximations of discontinuous functions. Yet, for instance, Gzyl and Palacios [23] have shown the absence of the Gibbs–Wilbraham phenomenon for Bernstein polynomials. It is still an open problem which properties of the approximation cause the Gibbs–Wilbraham phenomenon, but it is our conjecture that the order of accuracy for smooth functions is the deciding factor.

• •

While Bernstein polynomials converge uniformly for every continuous function, for instance, polynomial interpolation in general only converges if $u$ is at least (Dini–)Lipschitz continuous [9]. Even worse, Faber [16] showed in 1914 that no polynomial interpolation scheme – no matter how the points are distributed – will converge for the whole class of continuous functions. For approximations by orthogonal projection, such as Fourier series, a similar result follows from divergence of the corresponding operator norms [7, 44].

We conclude from Remark 3 that Bernstein polynomials, while appearing not attractive for approximating sufficiently smooth functions, provide some advantages for the approximation of just continuous or even discontinuous functions.

3 The Bernstein Procedure

In this section, we introduce a novel sub-cell shock capturing procedure by using Bernstein polynomials. The procedure is described for the reference element $[-1,1]$ and is based on replacing polluted high-order approximations by their Bernstein reconstruction.

3.1 Bernstein Reconstruction

We start by introducing the (modified) Bernstein reconstruction which is obtained by applying the Bernstein operator (6) to a polynomial approximation. Let $u\in\mathbb{P}_{N}([-1,1])$ be an approximate solution at time $t\geq 0$ with coefficients $\underline{\hat{u}}\in\mathbb{R}^{N+1}$ w. r. t. a (nodal or modal) basis $\{\varphi_{n}\}_{n=0}^{N}$ in the reference element $\Omega_{ref}=[-1,1]$. Then, the original approximation, for instance, obtained by interpolation or (pseudo) $L^{2}$-projection, is modified in the following way:

Compute the Bernstein reconstruction of $u$ by

[TABLE] 2. 2.

Write $B_{N}[u]$ w. r. t. the basis $\{\varphi_{n}\}_{n=0}^{N}$ by a change of bases with transformation matrix $\underline{\underline{T}}\,$, i. e.

[TABLE]

where $\underline{b}$ is a vector containing the Bernstein coefficients $\beta_{n,N}=u\left(-1+2\frac{n}{N}\right)$.

To put it in a nutshell, the idea is to replace the coefficients $\underline{\hat{u}}$ in a troubled cell by the coefficients $\underline{u}^{(B)}$ of the Bernstein reconstruction. Table 1 lists the condition numbers [56] of the transformation matrix $\underline{\underline{T}}\,$ for the Lagrange and Legendre bases. Thereby, the nodal basis of Lagrange polynomials is considered w. r. t. the Gauss–Lobatto points in $[-1,1]$. For all reasonable polynomial degrees ($N\leq 10$) and both bases we observe the condition number to be fairly small.

Concerning the now obtained Bernstein reconstruction, Lemma 1 does not only ensure preservation of bounds, but also allows us to enforce such bounds. Let us introduce the modified Bernstein reconstruction w. r. t. the lower bound $m$ and the upper bound $M$

[TABLE]

To the best of the authors’ knowledge, the modified Bernstein operator (29) has not been defined anywhere else yet. The most beautiful property of the modified Bernstein operator is preservation – by default – of lower and upper bounds $m$ and $M$, respectively. This follows directly from Lemma 1 and is summarised in

Theorem 4**.**

Let $B_{N}^{(m,M)}$ be the modified Bernstein operator w. r. t. the lower bound $m$ and the upper bound $M$ given by (29) and let $u$ be some function. Then $B_{N}^{(m,M)}[u]$ fulfils

[TABLE]

We close this subsection by noting that, for instance, positivity (think about density in the Euler equations), can be easily ensured by setting $m=\varepsilon$, where $\varepsilon>0$ is a suitable value larger than machine precision.

3.2 Proposed Procedure

Finally, we propose a procedure on how to replace the original polynomial approximation by its (modified) Bernstein reconstruction. Let $\Omega_{i}$ be a troubled element with approximate solution $u\in\mathbb{P}_{N}([-1,1])$ and let the coefficients of $u$ w. r. t. a (nodal or modal) basis $\{\varphi_{n}\}_{n=0}^{N}$ be given by $\underline{\hat{u}}\in\mathbb{R}^{N+1}$. The Bernstein procedure consists of two steps:

Compute the (modified) Bernstein reconstruction $B_{N}[u]$ of $u$. 2. 2.

Build an ’appropriate’ convex combination $u^{(\alpha)}$ of the original approximation $u$ and its Bernstein reconstruction $B_{N}[u]$, i. e.

[TABLE]

with $\alpha\in[0,1]$.

W. r. t. the original basis $\{\varphi_{n}\}_{n=0}^{N}$, and therefore utilising the transformation (28), the $\alpha$ Bernstein reconstruction $u^{(\alpha)}\in\mathbb{P}_{N}([-1,1])$ can be written as

[TABLE]

Hence, its coefficients are given by

[TABLE]

This procedure can be incorporated easily into an already existing solver. Note that we have $u^{(1)}=u$ and $u^{(0)}=u^{(B)}$ in the extreme cases of $\alpha=1$ and $\alpha=0$. Thus, the $\alpha$ Bernstein reconstruction $u^{(\alpha)}$ (linearly) varies between the original (and potentially oscillating or boundary violating) approximation $u$ and the more robust (modified) Bernstein reconstruction $B_{N}[u]$. This is illustrated in Figure 1 for the signum function $u(x)=\text{sign}(x)$, different parameters $\alpha$, and $N=1,5,9$.

Obviously, the order of the approximation is reduced to one in elements where $\alpha<1$. Yet, in elements with $\alpha=1$, which define the large majority of the domain, high accuracy is retained. Subsequently, we demonstrate how the parameter $\alpha$ can be adapted to the regularity of the underlying solution.

3.3 Selection of Parameter $\alpha$

The parameter $\alpha$ in (31) is adapted to adequately capture different discontinuities and regions of smoothness in the solution. The value of $\alpha$ adjusts in space and time to accurately capture strong variations in the solution. Hence, the proposed Bernstein procedure is able to calibrate the polynomial approximation to the regularity of the solution. It should be stressed that modification of the mesh topology, the number of degrees of freedom, node positions, or the type of SE method is utterly unnecessary. This makes the Bernstein procedure an efficient and easy to implement shock capturing method. As described above, the extreme values $\alpha=1$ and $\alpha=0$ yield the original approximation $u$ and the Bernstein reconstruction $B_{N}[u]$, respectively. For intermediate parameter values $\alpha\in(0,1)$, the Bernstein reconstruction renders the convex combination $u^{(\alpha)}=\alpha u+(1-\alpha)u^{(B)}$ more and more robust. Thus, $\alpha$ allows us to adapt the amount of stabilisation introduced by the Bernstein reconstruction. Here, we chose $\alpha$ being a function of a discontinuity sensor as proposed in [43, 27] and briefly revisited in the next subsection.

3.4 Discontinuity Sensor

In this work, we use a discontinuity sensor proposed in [17] to detect troubled elements and to steer the parameter $\alpha$. This sensor is an element-based function leading to a single scalar measure of the solutions’ smoothness. It is a nonlinear functional

[TABLE]

which depends on a sensing variable $s$. For systems, such as the Euler equations, the density or Mach number could be utilised for the sensing variable. Since we only consider scalar conservation laws in this work, $s$ is simply chosen as the conserved quantity $u$. The sensor was first proposed in [17] and is based on comparing polynomial annihilation (PA) operators [3] of increasing orders. A PA operator of order $m$,

[TABLE]

is designed as a high-order approximation of the jump function

[TABLE]

of the sensing variable $s$. Here, $S_{x}=\{\xi_{0}(x),\dots,\xi_{m}(x)\}\subset[-1,1]$ denotes a set of $m+1$ local grid points around $x$ and the so called annihilation coefficients $c_{j}\mathrel{\mathop{\ordinarycolon}}[-1,1]\to\mathbb{R}$ are given by

[TABLE]

where $\{p_{l}\}_{l=0}^{m}$ is a basis of $\mathbb{P}_{m}$. An explicit formula for the annihilation coefficients is provided in [3] by

[TABLE]

The normalisation factor $q_{m}$ is chosen as

[TABLE]

and ensures convergence to the right jump strengths. It has been proved in [3] that

[TABLE]

holds, where $h(\xi)=\max\{|\xi_{j}-\xi_{j-1}|\ |\ \xi_{j-1},\xi_{j}\in S_{x}\}$ and $I_{x}$ is the smallest closed interval such that $S_{x}\subset I_{x}$. We note that PA might be further enhanced, for instance, by applying a minmod limiter as described in [3].

Next, let the sensor value of order $m$ be given by

[TABLE]

i. e. by the greatest absolute value of the PA operator of order $m$ at the mid points of the collocation points. If $s$ has at least $m$ continuous derivatives, the PA operator $L_{m}$ provides convergence to [math] of order $m+1$, and we therefore expect the sensor value (41) to decrease for an increasing order $m$. In this case, a parameter value $\alpha=1$ ($\alpha>0$) is chosen, which means that the Bernstein procedure is not (fully) activated. In this work, we only compare the sensor values of order $m=1$ and $m=3$ since it was endorsed in [3] to use the same number of local grid points $\xi_{j}$ on both sides of a point $x$. Of course, a variety of modification is possible for the PA sensor. Now, the parameter value $\alpha=0$ is chosen only if

[TABLE]

holds, i. e. if the sensor value does not decrease when going over from order $m=1$ to $m=3$. Finally, we decide for the parameter $\alpha$ to linearly vary between $\alpha=1$ and $\alpha=0$. This is realised by the parameter function

[TABLE]

where $\kappa\in(0,1)$ is a problem dependent ramp parameter. Figure 2 illustrates the parameter function w. r. t. the values of the PA sensor $S$ given by (42).

For the later numerical tests we investigated different other parameter functions as well, of which some have been discussed in [27]. Yet, we obtained the best results with (43). It should be noted that the above revisited PA sensor is only recommended for high orders $N\geq 4$. In our numerical test, we have observed some miss-identifications for $N=3$. For $N\leq 3$ other discontinuity sensors could be used instead. We mention the modal-decay based sensor of Persson and Peraire [43, 27] and its refinements [4, 31] as well as the KXRCF sensor [32, 46] of Krivodonova et al., which is build up on a strong superconvergence phenomena of the DG method at outflow boundaries. Future work will address a detailed comparison of different shock sensors.

Remark 5**.**

The above sensor is simple to implement, but has its price, which is the introduction of the problem dependent and tunable parameter $\kappa\in(0,1)$ in (43). In practice, we follow a strategy proposed by Guermond, Pasquetti, and Popov for their entropy viscosity method for nonlinear conservation laws [22]. For a fixed polynomial degree $N$, the parameter $\kappa$ is tuned by testing the method on a coarse grid. For all problems presented later in §5, the tuning has been done quickly on a coarse mesh of $I=10$ elements. Further, we have observed the Bernstein procedure to be robust w. r. t. the parameter $\kappa$. This will be demonstrated in §5.1.1 for the linear advection equation.

4 Entropy, Total Variation, and Monotone (Shock) Profiles

In this section, we investigate some analytical properties of the proposed Bernstein procedure, such as its effect on entropy, total variation and monotone (shock) profiles of the numerical solution. For sake of simplicity, all subsequent investigations are carried out on the reference element $\Omega_{\mathrm{ref}}=[0,1]$.

4.1 Entropy Stability

Let $U\in C^{1}$ be a convex entropy function. We show that the proposed Bernstein reconstruction only yields a change in the total amount of entropy which is consistent with the total amount of the entropy of the original approximation $u$. This means that for increasing $N$ the total amount of entropy of the Bernstein reconstruction $B_{N}[u]$ converges to the total amount of entropy of the function $u$. Yet, even though the change of entropy is consistent, it does not always yield a decrease of the entropy. This is demonstrated by

Example 6**.**

Let $u(x)=x^{2}$ with Bernstein reconstruction $B_{N}[u](x)=x^{2}+\frac{1}{N}x(1-x)$. For the usual $L^{2}$-entropy $U(u)=u^{2}$, we have

[TABLE]

Thus, the total amount of entropy is increased by applying the Bernstein procedure.

In fact, Example 6 is no exceptional case. We note that for every continuous and convex function $u$, the sequence of Bernstein reconstructions will converge to $u$ from above, i. e.

[TABLE]

holds for all $N\geq 1$; see for instance [45, Theorem 7.1.8 and 7.1.9]. Hence, for $u\geq 0$, the $L^{2}$-entropy will be increased by the Bernstein procedure. Yet, we can further prove that the change of the total amount of entropy by the Bernstein reconstruction is consistent and vanishes for increasing $N$.

Theorem 7**.**

Let $U\in C^{1}$ be a convex entropy function and let $u$ be piecewise $C^{2}$. Then,

[TABLE]

holds.

Proof.

Since $U$ is $C^{1}$ and $u$ is piecewise $C^{2}$, equation (15) yields

[TABLE]

where $0<x_{1}<\dots<x_{k}<1$ are the points where $u^{\prime\prime}$ does not exist. Hence, $U\circ B_{N}[u]$ converges almost everywhere to $U\circ u$. Further, $u$ and all $B_{N}[u]$ are uniformly bounded, let us say by $m$ and $M$, i. e.

[TABLE]

for all $x\in[0,1]$. Since $U$ is continuous, $U$ is also bounded on $[m,M]$ and there is a $v^{*}\in[m,M]$ such that

[TABLE]

Thus, we have

[TABLE]

and the sequence $\left(U\circ B_{N}[u]\right)_{N\in\mathbb{N}}$ is uniformly bounded. Finally, Lebesgue’s dominated convergence theorem [15, Chapter 1.3] yields

[TABLE]

and therefore the assertion. ∎

Note that in the proposed Bernstein procedure, the Bernstein reconstruction is always computed for $u\in\mathbb{P}_{N}$. Yet, Theorem 7 holds for general $u$ which are piecewise $C^{2}$.

4.2 Total Variation

A fundamental property of the (exact) solution of a scalar hyperbolic conservation law (1), assuming the initial data function $u_{0}(x)=u(0,x)$ has bounded variation, is that [35, 55]

no additional spatial local extrema occur. 2. 2.

the values of local minima do not decrease and the values of local maxima do not increase.

As a consequence, the total variation (TV) of the solution,

[TABLE]

is a non-increasing function in time, i. e.

[TABLE]

for $t_{2}\geq t_{1}$. In the presence of (shock) discontinuities, in fact, the TV typically decreases [24]. We are thus interested in designing shock capturing methods which mimic this behaviour of being TV diminishing (TVD). The proposed Bernstein procedure is now shown to fulfil the TVD property in the sense that the Bernstein reconstruction $B_{N}[u]$ of a function $u$ has a reduced (or the same111In all numerical tests, we actually observed the TV to decrease) TV, i. e.

[TABLE]

holds for the Bernstein reconstruction (27).

Theorem 8**.**

Let $B_{N}[u]\in\mathbb{P}_{N}$ be the Bernstein reconstruction of a function $u\mathrel{\mathop{\ordinarycolon}}[0,1]\to\mathbb{R}$. Then, the TV of $B_{N}[u]$ is less or equal to the TV of $u$, and the Bernstein procedure is TVD.

Proof.

Since $B_{N}[u]\in\mathbb{P}_{N}$ and by consulting (14), we have

[TABLE]

where the integrals are given by

[TABLE]

Thus, inequality

[TABLE]

follows, and therefore the assertion. ∎

4.3 Monotone (Shock) Profiles

Let $u\mathrel{\mathop{\ordinarycolon}}[0,1]\to\mathbb{R}$ be a piecewise smooth function with single discontinuity at $x_{1}\in(0,1)$, representing a shock profile in a troubled element. It is desirable for the polynomial approximation of such a function to not introduce new (artificial) local extrema. Yet, typical polynomial approximations, such as interpolation and (pseudo) projections, are doing so by the Gibbs–Wilbraham phenomenon [26]. The Bernstein reconstruction, however, has been proved to not feature such spurious oscillations. This can, for instance, be noted from

Theorem 9** (Theorem 1.9.1 in Lorentz [39]).**

Suppose that $u$ is bounded in $[0,1]$ and let $L^{+},L^{-}$ respectively denote the right and left upper limits and $l^{+},l^{-}$ the right and left lower limits of $u$ at a point $x$. Then

[TABLE]

Note that Theorem 9 is fairly general. The absence of the Gibbs–Wilbraham phenomenon (without taking convergence into account) could have been noted by Lemma 1 already. It should be stressed that we can not just rule out spurious (Gibbs–Wilbraham) oscillations for the Bernstein reconstruction of discontinuous (shock) profiles, but we are further able to ensure the preservation of monotonicity. Let $u$ be a monotonic increasing (and possibly discontinuous) function on $[0,1]$. Then, the following lemma ensures that the Bernstein reconstruction $B_{N}[u]$ is monotonic increasing as well.

Lemma 10**.**

Let $u\mathrel{\mathop{\ordinarycolon}}[0,1]\to\mathbb{R}$ be monotonic increasing and let $B_{N}[u]$ denote the Bernstein reconstruction of $u$. Then, $B_{N}[u]$ is monotonic increasing as well.

Proof.

Note that for monotonic increasing $u$, we have

[TABLE]

Thus, by consulting (14), inequality

[TABLE]

follows and as a result also the assertion. ∎

Note that the same result holds for monotonic decreasing (shock) profiles.

5 Numerical Tests

In this section, we test the Bernstein procedure incorporated into a nodal collocation-type discontinuous Galerkin finite element method (DGFEM) as described in [25]. For time integration, we have used the explicit TVD/SSP-RK method of third order using three stages (SSPRK(3,3)) given in [21] by Gottlieb and Shu: Let $u^{n}$ be the solution at time $t^{n}$, then the solution $u^{n+1}$ at time $t^{n+1}$ is obtained by

[TABLE]

where $L(u)$ is a discretisation of the spatial operator of the DGFEM. For the time step size, we have used

[TABLE]

with $C=0.1$ and where $\max{|f^{\prime}(u)|}$ is calculated for all $u$ between $\min_{x\in\Omega}u_{0}(x)$ and $\max_{x\in\Omega}u_{0}(x)$ in all our numerical tests. Following [12] in parts, four different problems are investigated for which the exact solutions can be calculated. Note that we assume periodic boundary conditions (BCs) in all numerical tests. This restriction is utterly unnecessary for the proposed Bernstein procedure and is only made in order to compactly provide reference solutions, which refer to the exact entropy solutions. Further, for every problem the local Lax–Friedrichs (Rusanov) flux

[TABLE]

has been applied, where $\lambda_{\text{max}}=\max_{u\in[u_{-},u_{+}]}\left|f^{\prime}(u)\right|$ is a locally determined viscosity coefficient based on maximum characteristics speed [36, Chapter 12.5].

5.1 Linear Advection Equation

Let us consider the linear advection equation

[TABLE]

on $\Omega=[0,1]$ with periodic BCs and a discontinuous initial condition (IC)

[TABLE]

By the method of characteristics, the solution at time $t\geq 0$ is given by

[TABLE]

The linear advection equation is the simplest PDE that can feature discontinuous solutions. Thus, the (shock capturing) method can be observed in a well-understood setting, isolated from nonlinear effects. Yet, the linear advection equation provides a fairly challenging example. Similar to contact discontinuities in Euler’s equations, discontinuities are not self-steeping, i. e. once such a discontinuity is smeared by the method, it can not be recovered to its original sharp shape [33].

5.1.1 Parameter Study

We start by investigating the ramp parameter $\kappa\in(0,1)$, which goes into the PA sensor (43) and steers the Bernstein procedure. Figure 3 demonstrates the effect of the ramp parameter $\kappa\in(0,1)$ on the results produced by the Bernstein procedure for a linear advection equation with discontinuous IC. The IC has thereby been evolved over time until $t=1$. We observe the Bernstein procedure to be fairly robust w. r. t. the ramp parameter $\kappa$. On coarse meshes, such as $I=10$ and $I=20$, slight differences can be observed between different parameter values. Yet, these differences become less significant when the mesh is refined. As mentioned in Remark 5, the tuning of the ramp parameter has been done quickly on a coarse mesh of $I=10$ elements for all problems presented.

5.1.2 Comparison with Usual Filtering in DG Methods

Next, we investigate the approximation properties of a DGFEM enhanced with the Bernstein procedure for the linear advection equation with discontinuous IC. Figure 4 illustrates the results at time $t=1$. Further, we compare our results with the DGFEM without any filtering and with a usual filtering technique, where the DG solution $u\in\mathbb{P}_{N}(\Omega_{i})$ is replaced by its mean,

[TABLE]

in a troubled element $\Omega_{i}$. Troubled elements are detected by a critical value $S(u)\geq 1$, see (42). Note that the usual filtering is therefore expected to be applied in less elements than the Bernstein procedure, since the Bernstein procedure is already activated for $S(u)>\kappa$, see (43). Yet, we still observe the usual filtering to smear the numerical solution around discontinuities considerably. Over time, this smearing yields the numerical solution to nearly become constant. This can be observed in Figure 5, where the results are further evolved in time until $t=10$. At the same time, the results of the Bernstein procedure remain in their relatively sharp shape, even near discontinuities. Yet, no oscillations are observed, in contrast to the results produced by the DGFEM without any filtering. It should be stressed that this test case is especially challenging for the usual filtering by mean values, since the initial discontinuities have traveled trough the domain several times until $t=10$ is reached. The following tests might provide a fairer comparison. An extensive smearing of the numerical solution by the usual filtering by mean values, especially compared to the Bernstein procedure, is always observed to some extent, however.

Remark 11**.**

In Figure 4(g) (and Figure 9(d)) we note stronger oscillations for the DGFEM with mean value filtering than without any filtering. Typically, this is not excepted. A reason for this behaviour might be that the mean value filtering, when activated by the above sensor, is neither ensured to be TVD nor to preserve the relation between boundary values at element interfaces (while $u_{-}<u_{+}$ holds for the original approximation, mean value filtering in one or both elements connected by the interface might result in a reverse relation $\overline{u}_{-}>\overline{u}_{+}$). Such behaviour could be prevented by local projection limiters which are steered by modified minmod function; see [12, 11].

5.2 Invicid Burgers’ Equation

Let us now consider the nonlinear invicid Burgers’ equation

[TABLE]

on $\Omega=[0,1]$ with smooth IC

[TABLE]

and periodic BCs. For this problem a shock develops in the solution when the wave breaks at time

[TABLE]

In the subsequent numerical tests, we consider the solution at times $t=2$ and $t=3$. The reference solutions have been computed using characteristic tracing, solving the implicit equation $u(t,x)=u_{0}(x-tu)$ in smooth regions. The jump location, separating these regions, can be determined by the Rankine–Hugoniot condition. Figure 6 illustrates the results at time $t=2$, at which the shock waves starts to arise. As a consequence, only slight differences are observed at this state. In particular, we can observe that the Bernstein procedure does not smear the solution in smooth regions, even when steep gradients arise. A slight smearing can be observed for the usual filtering by mean values around the location of the arising shock at $x=0.5$. Figure 7 illustrates the results at time $t=3$, for which the shock has fully developed and has already traveled through the whole domain once. While we observe oscillations for the DGFEM without filtering and smearing for usual filtering by mean values, the Bernstein procedure still provides sharp profiles for the discontinuous solution.

5.3 A Concave Flux Function

Next, we investigate the conservation law

[TABLE]

with a concave flux function $f(u)=u(1-u)$ on $\Omega=[0,2]$, periodic BCs, and a discontinuous IC

[TABLE]

For this problem a rarefaction wave develops in the solution. Figure 8 illustrates the results at time $t=0.5$. Here, the DGFEM without any filtering fails to capture the rarefaction wave around $x=1$, yielding a wrong (physically unreasonable) weak solution. The usual filtering by mean values captures the rarefaction wave, but again smears the solution. The Bernstein procedure is also able to capture the rarefaction wave, while providing notably sharper profiles. Yet, we can observe some remaining slight oscillation for the Bernstein procedure in some cases. Remember that the Bernstein procedure is ensured to preserve bounds. Thus, the oscillations probably have been caused by the Bernstein procedure getting activated after oscillations have already been developed in the numerical solution by the DGFEM. Future work will investigate the possibility of using other shock sensors to steer the proposed Bernstein procedure in greater detail. These sensors might be activated once the TV increases over time or once the numerical solution starts to violate certain bounds.

5.4 The Buckley–Leverett Equation

Finally, we consider the Buckley–Leverett equation

[TABLE]

with a non-convex flux function $f(u)=\frac{u^{2}}{u^{2}+(1-u)^{2}}$ on $\Omega=[0,2]$, periodic BCs, and a discontinuous IC

[TABLE]

The Buckley–Leverett equation is often used to describe an immiscible displacement process, such as the displacement of oil by water [47, Chapter 4.2]. Due to its nonconvex flux function, the Riemann solution involves a so-called compound wave (sometimes also referred to as a composite curve), which contains a shock discontinuity and a rarefaction wave at the same time. For a nonlinear system of equations this can arise in any nonlinear field that fails to be genuinely nonlinear; see [55, Chapter 16.1] and [59, 60, 37]. The reference solution has been computed using characteristic tracing again. Figure 9 illustrates the results at time $t=0.5$. For this problem, the DGFEM without filtering provides fairly poor results, polluted by heavy oscillations, in all tests. In some cases for $N=3$ even the usual filtering by mean values shows some oscillations. This problem seems to be related to some miss-identifications of troubled elements by the PA sensor for $N=3$. For finer meshes the problem vanishes in case of the usual filtering by mean values. Yet, we still observe the Bernstein procedure to also perform relatively poor for $N=3$. For $N>3$, the Bernstein procedure is again observed to provide notably better results than usual filtering by mean values. We conclude that the Bernstein procedure steered by the PA sensor (43) can only be recommended for higher orders $N\geq 4$. Yet, similar observations for the usual filtering by mean values indicate that the relatively poor behaviour for $N=3$ might be caused by the PA sensor, rather than by the Bernstein procedure itself.

6 Concluding Thoughts

In this work, we have proposed a novel shock capturing procedure for SE approximations for scalar hyperbolic conservation laws. The procedure is easy to implement and neither increases the computational complexity of the method nor adds additional time step or CFL restrictions, as other common shock capturing methods do. The procedure essentially consists of going over from the original (oscillatory) polynomial approximation to its Bernstein reconstruction in troubled elements. Thereby, the Bernstein reconstruction of a function is obtained by applying the (modified) Bernstein operator (29). This operator has been proved to reduce the TV of the underlying function as well as to preserve monotone (shock) profiles. Both properties distinguish the shock capturing procedure as especially suitable for scalar conservation laws, which come along with a TVD property for their physically reasonable solutions. We further note that the modified procedure is not just able to preserve but even to enforce certain bounds for the numerical solution, such as positivity. Finally, the proposed procedure can be calibrated to the regularity of the underlying solution by using a discontinuity sensor which, in this work, is based on comparing the values of PA operators of increasing orders. Numerical tests demonstrate that the procedure is able to significantly enhance SE approximations in the presence of shocks.

Future work will focus on the extension of the proposed Bernstein procedure to unstructured meshes in multiple dimensions and systems of conservation laws. Further, the investigation of other shock sensors would be of great interest. These sensors might be activated once the TV increases over time or once the numerical solution starts to violate certain bounds.

Acknowledgements

The author would like to thank Thomas Sonar (TU Braunschweig, Germany) for helpful advice. This work is supported by the German Research Foundation (DFG, Deutsche Forschungsgemeinschaft) under Grant SO 363/15-1.