High Order Edge Sensors with $\ell^1$ Regularization for Enhanced Discontinuous Galerkin Methods

TL;DR

This paper introduces a novel high order edge sensor with $$ regularization for discontinuous Galerkin methods, improving the handling of troubled regions in hyperbolic conservation laws efficiently.

Contribution

It develops a high order edge sensor using polynomial annihilation and applies $$ regularization selectively to troubled cells in DG methods, enhancing accuracy and efficiency.

Findings

Effective detection of troubled elements with the edge sensor.

Improved solution accuracy in discontinuous regions.

Efficient implementation using ADMM for $$ optimization.

Abstract

This paper investigates the use of $\ell^1$ regularization for solving hyperbolic conservation laws based on high order discontinuous Galerkin (DG) approximations. We first use the polynomial annihilation method to construct a high order edge sensor which enables us to flag troubled elements. The DG approximation is enhanced in these troubled regions by activating $\ell^1$ regularization to promote sparsity in the corresponding jump function of the numerical solution. The resulting $\ell^1$ optimization problem is efficiently implemented using the alternating direction method of multipliers. By enacting $\ell^1$ regularization only in troubled cells, our method remains accurate and efficient, as no additional regularization or expensive iterative procedures are needed in smooth regions. We present results for the inviscid Burgers' equation as well as a nonlinear system of conservation laws using a nodal collocation-type DG method as a solver.

Tables

Table 1. Table 1 : Shock sensor for function displayed in Figure 2 .

Element	${\Omega }_1=(0,2)$	${\Omega }_2=(2,4)$	${\Omega }_3=(4,6)$	${\Omega }_4=(6,8)$	${\Omega }_5=(8,10)$
$S_1$	1.58	1.07	1.07	2.79	0.38
$S_3$	1.99	0.60	1.38	3.00	0.17
Discontinuity	yes	no	yes	yes	no

Table 2. Table 2 : Errors of the numerical solutions without and with ℓ 1 superscript ℓ 1 \ell^{1} regularization. In cases where the error value is NaN, the numerical solution broke down before the final time was reached. “no ℓ 1 superscript ℓ 1 \ell^{1} ” refers to the underlying DG method without ℓ 1 superscript ℓ 1 \ell^{1} regularization, “ ℓ 1 superscript ℓ 1 \ell^{1} ” refers to the DG method with ℓ 1 superscript ℓ 1 \ell^{1} regularization, and “ ℓ 1 superscript ℓ 1 \ell^{1} -mc” refers to the DG method with ℓ 1 superscript ℓ 1 \ell^{1} regularization and additional mass correction.

		${\parallel \cdot \parallel }_M$-error			${\parallel \cdot \parallel }_1$-error			${\parallel \cdot \parallel }_{\infty }$-error
$p$	$I$	no ${\ell }^1$	${\ell }^1$	${\ell }^1$-mc	no ${\ell }^1$	${\ell }^1$	${\ell }^1$-mc	no ${\ell }^1$	${\ell }^1$	${\ell }^1$-mc
3	15	1.5e-2	3.3e-2	3.3e-2	1.2e-2	2.2e-1	1.2e-2	1.7e-1	3.3e-1	3.3e-1
	31	1.5e-2	2.0e-2	2.0e-2	1.1e-2	1.2e-2	1.2e-2	2.2e-1	3.2e-1	3.2e-1
	63	1.9e-2	2.4e-2	2.4e-2	1.1e-2	1.1e-2 ✓	1.1e-2	3.6e-1	5.6e-1	5.6e-1
	127	3.3e-2	2.7e-2 ✓	2.7e-2	1.3e-2	1.1e-2 ✓	1.1e-2	1.4e-0	8.3e-1 ✓	8.3e-1
4	15	7.0e-2	5.9e-2 ✓	5.9e-2	3.9e-2	2.7e-2 ✓	2.7e-2	7.6e-1	7.5e-1 ✓	7.5e-1
	31	5.5e-2	4.6e-2 ✓	4.6e-2	2.4e-2	1.7e-2 ✓	1.7e-2	8.7e-1	8.6e-1 ✓	8.6e-1
	63	4.3e-2	3.9e-2 ✓	3.9e-2	1.6e-2	1.3e-2 ✓	1.3e-2	1.0	1.0e-0 ✓	1.0e-0
	127	3.8e-2	3.6e-2 ✓	3.6e-2	1.2e-2	1.2e-2 ✓	1.1e-2	1.4	1.3e-0 ✓	1.3e-0
5	15	1.6e-2	1.5e-2 ✓	1.5e-2	1.2e-2	1.2e-2 ✓	1.2e-2	2.5e-1	1.9e-1 ✓	1.7e-1
	31	2.0e-2	1.4e-2 ✓	1.3e-2	1.3e-2	1.0e-2 ✓	1.0e-2	5.3e-1	2.6e-1 ✓	2.5e-1
	63	7.3e-2	1.9e-2 ✓	1.6e-2	2.2e-2	1.0e-2 ✓	1.0e-2	3.6e-0	5.3e-1 ✓	4.3e-1
	127	NaN	2.7e-2 !	2.2e-2	NaN	1.0e-2 !	1.0e-2	NaN	1.1e-0 !	9.1e-1
6	15	5.9e-2	5.2e-2 ✓	5.2e-2	3.2e-2	2.3e-2	2.3e-2	8.1e-1	8.0e-1 ✓	8.0e-1
	31	4.7e-2	4.3e-2 ✓	4.3e-2	2.0e-2	1.6e-2 ✓	1.6e-2	9.6e-1	9.5e-1 ✓	9.5e-1
	63	3.9e-2	3.7e-2 ✓	3.6e-2	1.2e-2	1.2e-2 ✓	1.2e-2	1.2e-0	1.1e-0 ✓	1.1e-0
	127	NaN	3.4e-2 !	3.2e-2	NaN	1.1e-2 !	1.1e-2	NaN	1.3e-0 !	1.3e-0
7	15	1.7e-2	1.8e-2	1.8e-2	1.4e-2	1.2e-2 ✓	1.2e-2	2.2e-1	2.7e-1	2.8e-1
	31	2.5e-2	1.7e-2 ✓	1.8e-2	1.5e-2	1.1e-2 ✓	1.1e-2	7.0e-1	3.8e-1 ✓	4.0e-1
	63	5.0e-1	2.1e-2 ✓	2.2e-2	1.4e-1	1.0e-2 ✓	1.0e-2	1.1e+1	6.9e-1 ✓	7.5e-1
	127	2.8e-2	2.5e-2 ✓	2.5e-2	1.0e-2	1.0e-2 ✓	1.0e-2	1.3e-0	1.2e-0 ✓	1.2e-0
8	15	5.9e-2	4.9e-2 ✓	4.9e-2	3.2e-2	2.1e-2 ✓	2.1e-2	9.9e-1	8.6e-1 ✓	8.6e-1
	31	4.5e-2	4.0e-2 ✓	4.0e-2	1.9e-2	1.4e-2 ✓	1.4e-2	1.0e-0	1.0e-0 ✓	1.0e-0
	63	3.8e-2	3.6e-2 ✓	3.7e-2	1.2e-2	1.2e-2 ✓	1.2e-2	1.4e-0	1.3e-0 ✓	1.3e-0
	127	1.5e-1	3.2e-2 ✓	3.1e-2	3.1e-2	1.1e-2 ✓	1.1e-2	6.9e-0	1.4e-0 ✓	1.4e-0
9	15	2.0e-2	1.9e-2 ✓	1.8e-2	1.5e-2	1.3e-2 ✓	1.3e-2	4.5e-1	3.0e-1 ✓	2.9e-1
	31	7.8e-2	1.9e-2 ✓	1.8e-2	2.5e-2	1.1e-2 ✓	1.1e-2	4.5e-0	4.5e-1 ✓	4.1e-1
	63	NaN	2.3e-2 !	2.2e-2	NaN	1.0e-2 !	1.0e-2	NaN	9.0e-1 !	8.3e-1
	127	NaN	2.7e-2 !	2.8e-2	NaN	1.0e-2 !	1.0e-2	NaN	1.4e-0 !	1.4e-0

Table 3. Table 3 : u ​ ( 0 , x ) = sin ⁡ ( 2 ​ π ​ x ) 𝑢 0 𝑥 2 𝜋 𝑥 u(0,x)=\sin\left(2\pi x\right) . Errors of the numerical solutions without and with ℓ 1 superscript ℓ 1 \ell^{1} regularization. “no ℓ 1 superscript ℓ 1 \ell^{1} ” refers to the underlying DG method without ℓ 1 superscript ℓ 1 \ell^{1} regularization, “ ℓ 1 superscript ℓ 1 \ell^{1} ” refers to the DG method with ℓ 1 superscript ℓ 1 \ell^{1} regularization, and “ ℓ 1 superscript ℓ 1 \ell^{1} -mc” refers to the DG method with ℓ 1 superscript ℓ 1 \ell^{1} regularization and additional mass correction.

		${\parallel \cdot \parallel }_M$-error			${\parallel \cdot \parallel }_1$-error			${\parallel \cdot \parallel }_{\infty }$-error
$p$	$I$	no ${\ell }^1$	${\ell }^1$	${\ell }^1$-mc	no ${\ell }^1$	${\ell }^1$	${\ell }^1$-mc	no ${\ell }^1$	${\ell }^1$	${\ell }^1$-mc
3	2	1.2e-0	1.2e-0	1.2e-0	1.5e-0	1.6e-0	1.6e-0	9.5e-1	9.8e-1	9.8e-1
	4	1.3e-1	4.4e-1	4.4e-1	1.4e-1	5.8e-1	5.8e-1	1.3e-1	4.1e-1	4.1e-1
	8	6.3e-3	6.3e-3	6.3e-3	7.0e-3	7.0e-3	7.0e-3	1.2e-2	1.2e-2	1.2e-2
	16	3.8e-4	3.8e-4	3.8e-4	3.8e-4	3.8e-4	3.8e-4	9.9e-4	9.9e-4	9.9e-4
4	2	3.4e-1	9.2e-1	9.2e-1	4.2e-1	9.7e-1	9.6e-1	3.8e-1	9.0e-1	8.8e-1
	4	7.8e-3	7.8e-3	7.8e-3	1.0e-2	1.0e-2	1.0e-2	1.2e-2	1.2e-2	1.2e-2
	8	4.2e-4	4.2e-4	4.2e-4	4.4e-4	4.4e-4	4.4e-4	1.2e-3	1.2e-3	1.2e-3
	16	1.3e-5	1.3e-5	1.3e-5	1.3e-5	1.3e-5	1.3e-5	4.4e-5	4.4e-5	4.4e-5
5	2	8.0e-2	1.0e-0	1.0e-0	1.0e-1	1.3e-0	1.3e-0	1.3e-1	8.4e-1	7.8e-1
	4	2.1e-3	2.1e-3	2.1e-3	2.3e-3	2.3e-3	2.3e-3	5.3e-3	5.3e-3	5.3e-3
	8	2.8e-5	2.8e-5	2.8e-5	2.9e-5	2.9e-5	2.9e-5	7.7e-5	7.7e-5	7.7e-5
	16	1.2e-6	1.2e-6	1.2e-6	1.5e-6	1.5e-6	1.5e-6	1.6e-6	1.6e-6	1.6e-6
6	2	2.1e-2	9.9e-1	9.9e-1	2.5e-2	1.1e-0	1.1e-0	5.1e-2	9.9e-1	9.9e-1
	4	7.5e-5	7.5e-5	7.5e-5	8.0e-5	8.0e-5	8.0e-5	1.8e-4	1.8e-4	1.8e-4
	8	6.0e-6	6.0e-6	6.0e-6	7.6e-6	7.6e-6	7.6e-6	7.5e-6	7.5e-6	7.5e-6
	16	7.3e-7	7.3e-7	7.3e-7	9.4e-7	9.4e-7	9.4e-7	7.4e-7	7.4e-7	7.4e-7
7	2	2.0e-3	2.0e-3	2.0e-3	2.0e-3	2.0e-3	2.0e-3	6.4e-3	6.4e-3	6.4e-3
	4	3.9e-5	3.9e-5	3.9e-5	4.9e-5	4.9e-5	4.9e-5	6.8e-5	6.8e-5	6.8e-5
	8	3.9e-6	3.9e-6	3.9e-6	5.0e-6	5.0e-6	5.0e-6	4.1e-6	4.1e-6	4.1e-6
	16	4.9e-7	4.9e-7	4.9e-7	6.3e-7	6.3e-7	6.3e-7	4.9e-7	4.9e-7	4.9e-7

Equations

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

$${\partial }_{t}u+{\partial }_{x}f\left(u\right)=0,x\in \Omega ,$$

$$u\left(t,\xi \right)\approx u_p\left(t,\xi \right)=\sum _{k=0}^p{u}_k\left(t\right){\ell }_k\left(\xi \right),$$

$${\ell }_k\left(\xi \right)=\prod _{\begin{array}{c}i=0\\ i\ne k\end{array}}^p\frac{\xi -{\xi }_i}{{\xi }_k-{\xi }_i}$$

$$f\left(u\right)\approx f_p\left(t,\xi \right)=\sum _{k=0}^p{f}_k\left(t\right){\ell }_k\left(\xi \right),$$

$$\int_{-1}^{1}\dot{u}_{p}\ell_{i}\mathrm{d}\xi+\Bigl{.}\left(f^{\mathrm{num}}\ell_{i}\right)\Bigr{|}_{-1}^{1}-\int_{-1}^{1}f_{p}\ell_{i}^{\prime}\mathrm{d}\xi,\quad i=0,\dots,p.$$

$${\int }_{-1}^{1}g\left(\xi \right)\mathrm{d}\xi \approx {\int }_{-1}^1\sum _{k=0}^{p}g\left({\xi }_k\right){\ell }_k\left(\xi \right)\mathrm{d}\xi =\sum _{k=0}^p{\omega }_{k}g\left({\xi }_k\right).$$

$${⟨u,v⟩}_{\underline{\underline{M}}}=\sum _{k=0}^p{\omega }_{k}u\left({\xi }_k\right)v\left({\xi }_k\right)={\underline{u}}^T\underline{\underline{M}}\underline{v}$$

$$\underline{\underline{M}}=\mathrm{diag}\left([{\omega }_0,\dots ,{\omega }_p]\right)$$

$$\underline{u}=[u_0,\dots ,u_p{]}^T,\underline{v}=[v_0,\dots ,v_p{]}^T.$$

$$\left\langle{\dot{u}_{p},\,\ell_{i}}\right\rangle_{\underline{\underline{M}}\,}=\left\langle{f_{p},\,\ell_{i}^{\prime}}\right\rangle_{\underline{\underline{M}}\,}-\Bigl{.}\left(f^{\mathrm{num}}\ell_{i}\right)\Bigr{|}_{-1}^{1},\quad i=0,\dots,p.$$

$$\underline{\underline{M}}\dot{\underline{u}}=\underline{\underline{D}}{}^T\underline{\underline{M}}\underline{f}-\underline{\underline{R}}{}^T\underline{\underline{B}}{\underline{f}}^{\mathrm{num}},$$

$$\underline{\underline{D}}={\left({\ell }_i^{\prime }\left({\xi }_k\right)\right)}_{k,i=0}^p,\underline{\underline{R}}=\left(\begin{array}{cccc}0& \dots & 0& 1\\ 1& 0& \dots & 0\end{array}\right),\underline{\underline{B}}=\left(\begin{array}{cc}-1& 0\\ 0& 1\end{array}\right)$$

$$\underline{\underline{M}}\dot{\underline{u}}=-\underline{\underline{M}}\underline{\underline{D}}\underline{f}-\underline{\underline{R}}{}^T\underline{\underline{B}}\left({\underline{f}}^{\mathrm{num}}-\underline{\underline{R}}\underline{f}\right)$$

$$\underset{v\in {\mathbb{P}}_p\left({\Omega }_{\mathrm{ref}}\right)}{\mathrm{argmin}}H\left(v\right)\text{s.t.}\Vert v-u_p\Vert =0.$$

$$\underset{v\in {\mathbb{P}}_p\left({\Omega }_{\mathrm{ref}}\right)}{\mathrm{argmin}}H\left(v\right)\text{s.t.}\Vert v-u_p\Vert <\sigma$$

$$\underset{v\in {\mathbb{P}}_p\left({\Omega }_{\mathrm{ref}}\right)}{\mathrm{argmin}}\left({\Vert v-u_p\Vert }_2^2+\lambda H\left(v\right)\right)$$

$$H\left(v\right)={\Vert L_{m}v\Vert }_1,$$

$$[u]\left(\xi \right)=u\left({\xi }^{+}\right)-u\left({\xi }^{-}\right)$$

$$L_m[u]\left(\xi \right)=\frac{1}{q_m\left(\xi \right)}\sum _{x_j\in S_{\xi }}{c}_j\left(\xi \right)u\left(x_j\right),$$

$$S_{\xi }=\{x_0\left(\xi \right),\dots ,x_m\left(\xi \right)\}\subset {\Omega }_{\mathrm{ref}}$$

$$\sum _{x_j\in S_{\xi }}{c}_j\left(\xi \right)p_l\left(x_j\right)=p_l^{\left(m\right)}\left(\xi \right),j=0,\dots ,m,$$

$$c_j\left(\xi \right)=\frac{m!}{{\omega }_j\left(S_{\xi }\right)}\text{with}{\omega }_j\left(S_{\xi }\right)=\prod _{\begin{array}{c}{x}_i\in S_{\xi }\\ i\ne j\end{array}}\left(x_j-x_i\right)$$

$$q_m\left(\xi \right)=\sum _{x_j\in S_{\xi }^{+}}{c}_j\left(\xi \right),$$

$$L_{m}[u]\left(\xi\right)=\left\{\begin{array}[]{ccl}[u]\left(x\right)+\mathcal{O}\left(h\left(\xi\right)\right)&,&\text{if }x_{j-1}\leq\xi,x\leq x_{j},\\ \mathcal{O}\left(h^{\min\left(m,k\right)}\left(\xi\right)\right)&,&\text{if }u\in C^{k}\left(I_{\xi}\right),\end{array}\right.$$

$$u_p^{\mathrm{spar}}=\underset{v\in {\mathbb{P}}_p\left({\Omega }_{\mathrm{ref}}\right)}{\mathrm{argmin}}\left(\frac{1}{2}{\Vert v-u_p\Vert }_2^2+\lambda {\Vert L_{m}v\Vert }_1\right)$$

$$u_p^{\mathrm{spar}}=\underset{v\in {\mathbb{P}}_p\left({\Omega }_{\mathrm{ref}}\right)}{\mathrm{argmin}}\left({\Vert L_{m}v\Vert }_1+\frac{\mu }{2}{\Vert v-u_p\Vert }_2^2\right),$$

$$S_m=\underset{0\le k\le p-1}{\max }\left|L_m[u_p]\left({\xi }_{k+\frac{1}{2}}\right)\right|,$$

$${\xi }_{k+\frac{1}{2}}=\frac{{\xi }_k+{\xi }_{k+1}}{2},\text{for }k=0,\dots ,p-1,$$

$$S(u):=\frac{S_3}{S_1}\ge 1,$$

Full text

1 Introduction

This work is about a novel shock capturing procedure in spectral element (SE) type of methods for solving time-dependent hyperbolic conservation laws

[TABLE]

with smooth and discontinuous solutions. Such methods include the spectral difference (SD) [1, 2], discontinuous Galerkin (DG) [3], and flux reconstruction [4] methods. While all these methods provide fairly high orders of accuracy for smooth problems, they often lack desired stability and robustness properties, especially in the presence of (shock) discontinuities.

Many shock capturing techniques have therefore been developed over the last few decades. Such efforts date back more than $60$ years to the pioneering work of von Neumann and Richtmyer, [5] in which they add artificial viscosity terms to (1) in order to construct stable finite difference schemes for the equations of hydrodynamics. Since then, artificial viscosity has been added to a variety of algorithms [6, 7, 8, 9, 10, 11]. However, augmenting (1) with additional (higher) viscosity terms requires care about their design and size. Otherwise, new time stepping constraints for explicit methods can considerably decrease computational efficiency; see, e.g. [8, equation (2.1)]. Other interesting alternatives are based on (modal) filters [12, 13, 10] or applying viscosity to the different spectral scales [14, 15]. Finally, we mention those methods based on order reduction [16], mesh adaptation [17], and weighted essentially nonoscillatory (WENO) concepts [18, 19]. Yet, a number of issues regarding the fundamental convergence properties for these methods still remain unresolved. Moreover, even when the extension to multiple dimensions is straightforward, these schemes may be too computationally expensive for practical usage.

In this work we propose $\ell^{1}$ regularization as a novel tool to capture shocks in SE methods by promoting sparsity in the jump function of the approximate solution. $\ell^{1}$ regularization methods are frequently encountered in signal processing and imaging applications. They are still of limited use in solving partial differential equations numerically, however, and only a few studies, (see, e.g., [20, 21, 22, 23, 24, 25]) have considered sparsity or $\ell^{1}$ regularization of the numerical solution. A brief discussion of these investigations can be found in [24]. We note that while problems with discontinuous initial conditions were studied in [22, 23, 25], problems that form shocks were not. The technique developed in [20] was designed to promote sparsity in the frequency domain, making it less amenable to problems emitting shocks, where the frequency domain is not sparse. Further, with the exception of [24], each of these investigations applied $\ell^{1}$ regularization directly to numerical solution or to its residual, rather than incorporating it directly in the time stepping evolution.

In this investigation we follow the approach in [24], which incorporates $\ell^{1}$ regularization directly into the time dependent solver. Specifically, we promote the sparsity of the jump function that corresponds to the discontinuous solution. The jump function approximation is performed using the (high order) polynomial annihilation (PA) operator [26, 27, 28], which eliminates the unwanted “staircasing” effect, a common degradation of detail of the piecewise smooth solution arising from the classical total variation (TV) regularization. More specifically, the high order PA operator allows the resulting solution to be comprised of piecewise polynomials instead of piecewise constants. We solve the resulting $\ell^{1}$ optimization problem by the alternating direction method of multipliers (ADMM) [29, 30, 31]. A similar application of $\ell_{1}$ regularization was used in [24] to numerically solve hyperbolic conservation laws, though only for the Lax–Wendroff scheme and Chebyshev and Fourier spectral methods. It was concluded in [24] that although the Lax–Wendroff scheme yielded sufficient accuracy for relatively simple problems, its lower order convergence properties made it difficult to resolve more complicated ones. The new technique fared better using Chebyshev polynomials, although their global construction made it difficult to resolve the local structures without oscillations or excessive smoothing.

One possible solution is to use SE methods as the underlying mechanism for solving the hyperbolic conservation law. SE methods have the advantage of being more localized, for instance, and allow element-to-element variations in the optimization problem. In particular for the method we develop here, $\ell^{1}$ regularization is only activated in troubled elements, i.e., in elements where discontinuities are detected. This further enhances efficiency of the method. In the process, a novel discontinuity sensor based on PA operators of increasing orders is proposed, which is able to flag troubled elements. The discontinuity sensor steers the optimization problem and thus locally calibrates the method with respect to the smoothness of the solution. Numerical tests are performed for a nodal collocation-type DG method and the inviscid Burgers’ equation as well as for a nonlinear system. It should be stressed that the proposed procedure also carries over to other classes of methods, with the obvious extension to SE type methods. The extension to other types of methods, such as finite volume methods, is also possible under slight modifications of the procedure. We should also stress that in our new development of the $\ell^{1}$ regularization method for solving PDEs that admit shocks in their solutions that we do not require different methods to be used in smooth and nonsmooth regions. Such methods have been developed in, for instance, [32] and have been shown to be effective. Here we demonstrate that it is possible to avoid such additional complexities.

The rest of this paper is organized as follows: §2 briefly reviews the nodal collocation-type DG method, $\ell^{1}$ regularization, and PA operators which are needed for the development of our method. In §3 we describe the application of $\ell^{1}$ regularization by higher-order edge detectors to SE type methods. Further, a novel discontinuity sensor based on PA operators of increasing orders is proposed. Numerical tests for the inviscid Burgers’ equation, the linear advection equation, and a nonlinear system of conservation laws are presented in §4. The tests demonstrate that we are able to better resolve numerical solutions when $\ell^{1}$ regularization is utilized. We close this work with concluding thoughts in §5.

2 Preliminaries

In this section we briefly review all necessary concepts in order to introduce $\ell^{1}$ regularisation into the framework of discontinuous Galerkin methods in the subsequent section.

2.1 A nodal discontinuous Galerkin method

Let us consider a hyperbolic conservation law

[TABLE]

with suitable initial and boundary conditions. The domain $\Omega\subset\mathbb{R}$ is decomposed into $I$ disjoint, face-conforming elements $\Omega_{i}$, $\Omega=\bigcup_{i=1}^{I}\Omega_{i}$. All elements are mapped to a reference element, typically $\Omega_{\mathrm{ref}}=[-1,1]$, where all computations are performed.

In this work we consider a nodal collocation-type discontinuous Galerkin (DG) method [3] on the reference element. The solution $u$ as well as the flux function $f$ are approximated by interpolation polynomials of the same degree, giving the advantage of highly efficient operators. We further collocate the flux approximation based on interpolation with the numerical quadrature used for the evaluation of the inner products [33].

The first step is to introduce a nodal polynomial approximation

[TABLE]

where $p$ is the polynomial degree and $\{u_{k}\}_{k=0}^{p}$ are the $p+1$ time dependent nodal degrees of freedom at the element grid nodes $-1\leq\xi_{0}<\dots<\xi_{p}\leq 1$. Common choices are either the Gauss-Legendre or Gauss-Lobatto nodes [34]. We use Gauss-Lobatto points in the latter numerical tests, as they include the boundary points $\xi_{0}=-1,\xi_{p}=1$ and thus render the method more robust [35, 36, 37]. For the procedure of $\ell^{1}$ regularisation proposed in this work the specific choice of grid nodes is not crucial however.

Further, the Lagrange basis functions of degree $p$ are given by

[TABLE]

and satisfy the cardinal property $\ell_{k}\left(\xi_{i}\right)=\delta_{ki}$. The flux function is approximated in the same way, i.e.

[TABLE]

where, collocating the nodes for both approximations, the nodal degrees of freedom $\{f_{k}\}_{k=0}^{p}$ are given by $f_{k}\left(t\right)=f\left(u_{k}\left(t\right)\right)$.

We now obtain the formulation of the nodal DG method by inserting the polynomial approximations (3) and (5) into the conservation law (2), multiplying by a test function $\ell\in\{\ell_{k}\}_{k=0}^{p}$, integrating over the reference element, and applying integration by parts, resulting in

[TABLE]

Here, $f^{\mathrm{num}}$ is a suitably chosen numerical flux, providing a mechanism to couple the solutions across elements [38]. Further, $\dot{u}_{p}$ denotes the time derivative of the approximation while $\ell_{i}^{\prime}$ denotes the spatial derivative of the basis element with respect to $\xi$.

Next, the integrals are approximated by an (interpolatory) quadrature rule using the same nodes,

[TABLE]

From this quadrature rule, we introduce a discrete inner product

[TABLE]

with mass matrix

[TABLE]

and vectors of nodal degrees of freedom

[TABLE]

Using the discrete inner product, the spacial approximation (6) becomes

[TABLE]

Finally going over to a matrix vector representation and utilising the cardinal property of the Lagrange basis functions, the DG approximation can be compactly rewritten in its weak form as

[TABLE]

where

[TABLE]

are the differentiation matrix, restriction matrix, and boundary matrix, respectively. Let $\underline{f}^{\mathrm{num}}$ denote the vector containing the values of the numerical flux at the element boundaries. Note that applying integration by parts a second time to (6) would result in the strong form

[TABLE]

of the DG approximations. Both forms are equivalent when using summation by parts operators, which satisfy ${\underline{\underline{M}}\,\underline{\underline{D}}\,+\underline{\underline{D}}\,^{T}\underline{\underline{M}}\,=\underline{\underline{R}}\,^{T}\underline{\underline{B}}\,\underline{\underline{R}}\,}$. The strong form (14) can further be recovered as a special case of flux reconstruction schemes [39, 40, 41].

2.2 $\ell^{1}$ regularisation

Let $u\left(\xi\right)=u\left(t,\xi\right)$ be the unknown solution on an element $\Omega_{i}$ transformed into the reference element $\Omega_{\mathrm{ref}}$ and $u_{p}\in\mathbb{P}_{p}\left(\Omega_{\mathrm{ref}}\right)$ a spacial polynomial approximation at fixed time $t$. Assume that some measurable features of $u$ have sparse representation. Consequently, the approximation $u_{p}$ is desired to have this sparse representation as well.

Let $H$ be a regularisation functional which measures sparsity. The objective is to then solve the constrained optimisation problem

[TABLE]

The equality constraint, referred to as the data fidelity term, measures how accurately the reconstructed approximation fits the given data with respect to some seminorm $\mathinner{\!\left\lVert\cdot\right\rVert}$. Typically, the continuous $L^{2}$-norm $\mathinner{\!\left\lVert f\right\rVert}_{2}^{2}=\int|f|^{2}$ or some discrete counterpart is used. The regularisation term $H\left(v\right)$ enforces the known sparsity present in the underlying solution $u$ by penalising missing sparsity in the approximation. The regularisation functional $H$ further restricts the approximation space to a desired class of functions, here $\mathbb{P}_{p}\left(\Omega_{\mathrm{ref}}\right)$. Note that any $p$-norm with $p\leq 1$ will enforce sparsity in the approximation. In this work, we choose $H$ to be the $\ell^{1}$-norm of certain transformations of $v$.

It should be stressed that if $\mathinner{\!\left\lVert\cdot\right\rVert}$ is not just a seminorm but a strictly convex norm, for instance induced by an inner product, the equality constraint immediately and uniquely determines the approximation. Thus, instead of (15), typically the related denoising problem

[TABLE]

with $\sigma>0$ is solved, which relaxes the equality constraint on the data fidelity term. Equivalently, the denoising problem (16) can also be formulated as the unconstrained (or penalised) problem

[TABLE]

by introducing a non-negative regularisation parameter $\lambda\geq 0$. $\lambda$ represents the trade-off between fidelity to the original approximation and sparsity.

The unconstrained problem (17) is often solved with $H\left(v\right)=TV\left(v\right)$, where $TV(v)$ is the total variation of $v$. Following [24] however, in this work we solve (17) with

[TABLE]

where $L_{m}$ is a polynomial annihilation operator introduced in the next subsection. Using higher order polynomial annihilation operators will help to eliminate the staircase effect that occurs when using the TV operator (polynomial annihilation for $m=1$) for $H$, see [42].

2.3 Polynomial annihilation

Polynomial annihilation (PA) operators were originally proposed in [27]. One main advantage in using PA operators of higher orders ($m>1$) for the regularisation functional $H$ is that they allow distinction between jump discontinuities and steep gradients, which is critical in the numerical treatment of nonlinear conservation laws. PA regularisation is also preferable to TV regularisation when the resolution is poor, even when the underlying solution is piecewise constant.

Let $u\left(\xi^{-}\right)$ and $u\left(\xi^{+}\right)$ respectively denote the left and right hand side limits of $u\mathrel{\mathop{\ordinarycolon}}\Omega_{\mathrm{ref}}=[a,b]\to\mathbb{R}$ at $\xi$. We define the jump function of $u$ as

[TABLE]

and note that $[u]\left(\xi\right)=0$ at every $\xi\in\Omega_{\mathrm{ref}}$ where $u$ has no jump. We thus say that the jump function $[u]$ has a sparse representation. The polynomial annihilation operator of order $m$,

[TABLE]

is designed in order to approximate the jump function $\left[u\right]$. Here

[TABLE]

is a set of $m+1$ local grid points around $\xi$, the annihilation coefficients $c_{j}\mathrel{\mathop{\ordinarycolon}}\Omega_{\mathrm{ref}}\to\mathbb{R}$ are given by

[TABLE]

and $\{p_{l}\}_{l=0}^{m}$ is a basis of $\mathbb{P}_{m}\left(\Omega_{\mathrm{ref}}\right)$. An explicit formula for the annihilation coefficients utilising Newton’s divided differences is given by ([27]):

[TABLE]

for $j=0,\dots,m$. Finally, the normalisation factor $q_{m}$, calculated as

[TABLE]

ensures convergence to the right jump strength at every discontinuity. Here $S_{\xi}^{+}$ denotes the set $\{x_{j}\in S_{\xi}\ |\ x_{j}\geq\xi\}$ of all local grid points to the right of $\xi$.

In this work, the PA operator is applied to the reference element $\Omega_{\mathrm{ref}}=[-1,1]$ of the underlying nodal DG method using $p+1$ collocation points $\{\xi_{k}\}_{k=0}^{p}$, typically Gauss-Lobatto points including the boundaries. We can thus construct polynomial annihilation operators up to order $p$ by allowing the sets of local grid points $S_{\xi}$ to be certain subsets of the $p+1$ collocation points.

In [27] it was shown that

[TABLE]

where $h\left(\xi\right)=\max\{|x_{i}-x_{i-1}\ |\ x_{i-1},x_{i}\in S_{\xi}\}$ and $I_{\xi}$ is the smallest closed interval such that $S_{\xi}\subset I_{\xi}$. Note that $h(\xi)$ depends on the density of set of local grid points around $\xi$. Thus, $L_{m}$ provides $m$th order convergence in regions where $u\in C^{m}$, and yields a first order approximation to the jump. It should be stressed that oscillations develop around points of discontinuity as $m$ increases. The impact of the oscillation could be reduced by post-processing methods, such as the minmod limiter, as was done in [27]. However, as long as there is sufficient resolution between two shock locations, such oscillations do not directly impact our method. This is because we use the PA operator not to detect the precise location of jump discontinuities, but rather to enforce sparsity.

Figures 1(b) and 1(c) demonstrate the PA operator for the discontinuous function illustrated in Figure 1(a), while Figures 1(e) and 1(f) do so for the continuous but not differentiable function shown in Figure 1(g), and Figures 1(h) and 1(i) illustrate the PA operator for the smooth function in Figure 1(g).

For each function we used eight Gauss-Lobatto points to compute the PA operators for $m=1,3$ at the mid points $\xi_{j+1/2}=\left(\xi_{j-1}+\xi_{j}\right)/2$. As illustrated in Figure 1(b) and 1(c), the jump function approximation becomes oscillatory as $m$ is increased from $1$ to $3$, especially in the region of the discontinuity. The maximal absolute value of the PA operator is also not decreasing in this case. On the other hand, for the continuous and smooth functions displayed in Figures 1(d) and 1(g), the maximal absolute value decreases significantly when going from $m=1$ to $m=3$. This is consistent with the results in (25), and will be used to construct the discontinuity sensor in §3.3. For a discussion on the convergence of the PA operator see [27].

Remark 1**.**

In this work, we only consider one dimensional conservation laws. It should be stressed, however, that polynomial annihilation can be extended to multivariate irregular data in any domain. It was demonstrated in [27] that polynomial annihilation is numerically cost efficient and entirely independent of any specific shape or complexity of boundaries. In particular, in [43] and [44] the method was applied to high dimensional functions that arise when solving stochastic partial differential equations, which reside in a high dimensional space which includes the original space and time domains as well as additional random dimensions.

3 Application of $\ell_{1}$ regularization

In this section, we describe how the proposed $\ell^{1}$ regularization using PA operators, i.e., $H\left(v\right)=\mathinner{\!\left\lVert L_{m}v\right\rVert}_{1}$ in (17), can be incorporated into a DG method. While this kind of regularization functional was already investigated in [24], this work is the first to extend these ideas to an SE method and thus to allow element-to-element variations in the optimization problem. It should be stressed that the subsequent procedure relies on a piecewise polynomial approximation in space. Yet, by appropriate modifications of the procedure, it is also possible to apply $\ell^{1}$ regularization (with PA) to any other method.

3.1 Procedure

One of the main challenges in solving nonlinear conservation laws (2) is balancing high resolution properties and the amount of viscosity introduced to maintain stability, especially near shocks [5, 7, 11]. Applying the techniques presented in §2.2 and §2.3, we are now able to adapt the nodal DG method described in §2.1 to include $\ell^{1}$ regularization.

Our procedure consists of replacing the usual polynomial approximation $u_{p}$ by a sparse reconstruction

[TABLE]

with regularization parameter $\lambda$ in troubled elements after every (or every $k$th) time step by an explicit time integrator.

For the ADMM described in §3.5, it is advantageous to rewrite (26) in the usual form of an $\ell^{1}$ regularized problem as

[TABLE]

where $\mu=\frac{2}{\lambda}$ is referred to as the data fidelity parameter. Note that (26) and (27) are equivalent. In the later numerical tests, the data fidelity parameter $\mu$ and the regularization parameter $\lambda$ will be steered by a discontinuity sensor proposed in §3.3.

Remark 2**.**

One of the main drawbacks in using $\ell^{1}$ regularization for solving numerical partial differential equations, as well as for image restoration or sparse signal recovery, is in choosing the regularization parameter $\lambda$ (or $\mu$). Ideally, one would want to balance the terms in (26) or (27), but this is difficult to do without knowing their comparative size. Indeed, the $\ell^{1}$ regularization term $\mathinner{\!\left\lVert L_{m}v\right\rVert}_{1}$ heavily depends on the magnitudes of nonzero values in the sparsity domain, in this case the jumps. Larger jumps are penalized significantly more in the $\ell^{1}$ norm than smaller values.

Iterative spacially varying weighted $\ell^{1}$ regularization techniques (see, e.g., [45, 46]) are designed to help reduce the size of the norm, since the remaining values should be close to zero in magnitude. Specifically, the jump discontinuities which are meant to be in the solution can “pass through” the minimization. In this way, with some underlying assumptions made on the accuracy of the fidelity term, one could argue that both terms are close to zero. Consequently, the choice of $\lambda$ (or $\mu$) should not have as much impact on the results, leading to greater robustness overall. For the numerical experiments in this investigation, we simply chose regularization parameters that worked well. We did not attempt to optimize our results and leave parameter selection to future work.

3.2 Selection of the regularization parameter $\lambda$

The $\ell^{1}$ regularization should only be activated in troubled elements. In particular, we do not want to unnecessarily degrade the accuracy in the smooth regions of the solution. We thus propose to adapt the regularization parameter $\lambda$ in (26) to appropriately capture different discontinuities and regions of smoothness. As a result, the optimization problem will be able to calibrate the resulting sparse reconstruction to the smoothness of the solution. More specifically, to avoid unnecessary regularization, we choose $\lambda=0$ in elements corresponding to smooth regions. Note that this also renders the proposed method more efficient.

On the other hand, when a discontinuity is detected in an element, $\ell^{1}$ regularization will be fully activated by choosing $\lambda=\lambda_{\text{max}}$ in (26), which corresponds to the amount of regularization necessary to reconstruct sharp shock profiles. While no effort was made to optimize or even adapt this parameter, we found that using $\lambda_{\text{max}}=4\cdot 10^{2}$ in all of our numerical experiments yielded good results. A heuristic explanation for choosing $\lambda_{\text{max}}$ in this way stems from the goal of balancing the size of $\mathinner{\!\left\lVert L_{m}v\right\rVert}_{1}$ with the expected size of the fidelity term, which in this case means to be consistent with the order of accuracy of the underlying numerical PDE solver. As mentioned previously, choosing an appropriate $\lambda$ will be the subject of future work.

Between these extreme cases, i.e., $\lambda=0$ and $\lambda=\lambda_{\text{max}}$, we allow the regularization parameter to linearly vary and choose $\lambda$ as a function of the discontinuity sensor proposed in §3.3. As a consequence, we obtain more accurate sparse reconstructions while still maintaining stability in regions around discontinuities.

3.3 Discontinuity sensor

We now describe the discontinuity sensor which is used to activate the $\ell^{1}$ regularization and to steer the regularization parameter $\lambda$ in (26). The sensor is based on comparing PA operators of increasing orders. To the best of our knowledge this is the first time the PA operator is utilized for shock (discontinuity) detection in a PDE solver.111Of course (W)ENO schemes [47, 18, 19] compare slope magnitudes for determining troubled elements and choosing approximation stencils, so in this regard our method was inspired by WENO type methods.

At least for smooth solutions, discontinuous Galerkin methods are capable of spectral orders of accuracy. $\ell^{1}$ regularization as well as any other shock capturing procedure [5, 7, 13, 10, 11] should thus be just applied in (and near) elements where discontinuities are present. We refer to those elements as troubled elements.

Many shock and discontinuity sensors have been proposed over the last $20$ years for the selective application of shock capturing methods. Some of them use information about the $L^{2}$-norm of the residual of the variational form [48, 49], the primary orientation of the discontinuity [50], the magnitude of the facial interelement jumps [51, 52], or entropy pairs [53] to detect troubled elements. Others are not just able to detect troubled elements, but even the location of discontinuities in the corresponding element, such as the concentration method in [54, 55, 56, 57, 58, 32, 59]. The PA operator is also capable of detecting the location of a discontinuity [27, 60, 43] up to the resolution size, as in (25), using nonlinear postprocessing techniques. It should be stressed, however, that in this work we do not fully make use of this feature.

In what follows we present a novel discontinuity sensor based on the PA operator. Let the sensor value of order $m$ be

[TABLE]

i.e., the greatest absolute value of the PA operator at the midpoints

[TABLE]

of the collocation points. Since $L_{m}$ provides convergence to [math] of order $m$ in elements where $u$ has $m$ continuous derivatives, we expect $S_{3}<S_{1}$ to hold for an at least continuous function. Thus, $\ell^{1}$ regularization should just be fully activated if

[TABLE]

i.e., if the sensor value does not decrease as $m$ increases from $1$ to $3$. Sensor values of order $1$ and $3$ (rather than $2$) are chosen since having symmetry of the grid points in $S_{\xi}$ surrounding the point $\xi$ yields a simpler form for implementation [27]. Various modifications of the PA sensor (30) are possible and will be the topic of future research. In the following, the resulting PA sensor is demonstrated for the function displayed in Figure 2 on the interval $[0,10]$. We decompose the interval into $I=5$ elements and apply the PA operator and resulting PA sensor separately on each element.

Table 1 lists the sensor values of order $m=1$ and $m=3$ on each element in Figure 2. The last row further shows if the PA sensor (30) reacts. As can be seen in Table 1, the discontinuity sensor exactly identifies the troubled elements, simply by comparing the sensor values of order $m=1$ and $m=3$.

Even though we use this sensor only to detect troubled elements, in principle, this new sensor can also be used to more precisely determine the location and strength of the discontinuities. Such information can be used for refined domain decomposition [43]. Finally, we decide for the regularization parameter $\lambda$ to linearly vary between $\lambda=0$ and $\lambda=\lambda_{\text{max}}$ and thus utilize the parameter function

[TABLE]

where $\kappa\in[0,1)$ is a problem dependent ramp parameter. Observe that using $\lambda(S)$ is comparable to employing the weighted $\ell_{1}$ regularization as discussed in Remark 2.

For the later numerical tests we also considered other parameter functions, some as discussed in [61]. Yet the best results were obtained with (31). The same holds for other discontinuity sensors, such as the modal-decay based sensor of Persson and Peraire [7, 61] and its refinements [51, 8] as well as the KXRCF sensor [62, 47] of Krivodonova et al., which is built up on a strong superconvergence phenomenon of the DG method at outflow boundaries. For brevity, those results are omitted here.

Remark 3**.**

We note that the PA sensor might produce false positive or false negative misidentifications in certain cases. A false negative misidentification might arise from a discontinuity where the solution is detected to be smooth. This is encountered by the ramp parameter $\kappa$, which is observed to work robustly for $\kappa=0.8$ or $\kappa=0.9$ in all later numerical tests. A false positive misidentification might arise from a smooth solution which is detected to be nonsmooth (possibly discontinuous). In this case smooth parts of the solution with steep gradients will result in significantly greater values of $L_{m}u$ than parts of the solution with less steep gradients. As a result, the standard $\ell^{1}$ regularization would heavily penalize these features of the solution, yielding inappropriate smearing of steep gradients in smooth regions. By making $\lambda=\lambda(S)$ dependent on the sensor value, this problem can be somewhat alleviated. Using a weighted $\ell^{1}$ regularization, as suggested in Remark 2, should also reduce the unwanted smearing effect. Failure to detect a discontinuity would, after a number of time steps, yield instability. However it is unlikely that this would occur as the growing oscillations would more likely be identified as shock discontinuities.

3.4 Efficient implementation of the PA operator

While the PA operator was defined on the interior of the reference element $\Omega_{\mathrm{ref}}=[-1,1]$ in §2.3, the shock sensor proposed in §3.3 only relies on the values of the PA operator at the $p$ midpoints $\{\xi_{k+\frac{1}{2}}\}_{k=0}^{p-1}$ of the collocation points $\{\xi_{k}\}_{k=0}^{p}$. The same holds for the $\ell^{1}$ regularization term $H\left(v\right)=\mathinner{\!\left\lVert L_{m}v\right\rVert}_{1}$. The $\ell^{1}$-norm of the PA transformation is thus given by

[TABLE]

where the vector $\underline{v}$ once more consists of nodal degrees of freedom. We now aim to provide an efficient implementation of the PA operator $L_{m}$ in form of a matrix representation $\underline{\underline{L_{m}}}\,$, which maps the nodal values $\underline{v}$ to the values of the PA operator at the midpoints. Revisiting (20), this matrix representation is given by

[TABLE]

with

[TABLE]

and

[TABLE]

where

[TABLE]

Utilizing all prior matrix vector representations, we can now give the discretization of the $\ell^{1}$ regularized optimization problem (27) by

[TABLE]

Thus, we are able to solve the optimization problem directly for the nodal degrees of freedom of the sparse reconstruction $u_{p}^{\mathrm{spar}}$. Alternatively, the fidelity term $\mathinner{\!\left\lVert v-u_{p}\right\rVert}$ can also be approximated as

[TABLE]

i.e., by the Euclidean norm, yielding

[TABLE]

instead of (41). Future works will investigate the influence of the choice of the discrete norm on the performance of the $\ell^{1}$ regularization. Here we decided to use the Euclidean norm and thus the minimization problem (43), making the computations in §3.5 more intelligible.

3.5 The alternating direction method of multipliers

Many techniques have been recently proposed to solve optimization problems in the form of (41). Following [24], we use the ADMM [29, 30, 31] in our implementation. The ADMM has its roots in [63] and details of its convergence properties can be found in [64, 65, 66]. In the context of $\ell^{1}$ regularization, ADMM is commonly implemented using the split Bregman method [67], which is known to be an efficient solver for a broad class of optimization problems. To implement the ADMM, it is first necessary to eliminate all nonlinear terms within the $\ell^{1}$-norm. We thus introduce a slack variable

[TABLE]

and formulate (41) equivalently as

[TABLE]

To solve (45), we further introduce Lagrangian multipliers $\underline{\sigma}\in\mathbb{R}^{p},\underline{\delta}\in\mathbb{R}^{p+1}$ and solve the unconstrained minimization problem given by

[TABLE]

with objective function

[TABLE]

Here, $\beta>0$ is an additional positive regularization parameter and recall that the data fidelity parameter $\mu$ is given by $\mu=\frac{2}{\lambda}$ for $\lambda>0$; see (26) and (27). Note that if the Lagragian multipliers $\underline{\sigma},\underline{\delta}$ are updated a sufficient number of times, the solution of the unconstrained problem (46) will converge to the solution of the constrained problem (45). In the ADMM, the solution is approximated by alternating between minimizations of $\underline{v}$ and $\underline{g}$. A crucial advantage of this method is that, given the current value of $\underline{v}$ as well as the Lagrangian multipliers, the optimal value of g can be exactly determined by the shrinkage-like formula [67]

[TABLE]

where

[TABLE]

Given the current value $\underline{g}_{k+1}$, on the other hand, the optimal value of $\underline{v}$ is computed by the gradient descent method as

[TABLE]

where

[TABLE]

and the step size $\alpha>0$ is chosen to provide a sufficient descent in direction of the gradient. Finally, the Lagrangian multipliers are updated after each iteration by

[TABLE]

Algorithm 1 is borrowed from [31, 24] and compactly describes the above ADMM.

Figure 3 demonstrates the effect of $\ell^{1}$ regularization to a polynomial approximation with degree $p=13$ of a discontinuous sawtooth function $u(x)=\operatorname{sign}(x)-x$ on $[-1,1]$.

Due to the Gibbs phenomenon, the original polynomial approximation $u_{p}$ in Figure 3(a) shows spurious oscillations [68]. As a consequence, the approximation of the corresponding jump function by the PA operator for $m=3$ in Figure 3(b) also shows undesired oscillations away from the detected discontinuity at $x=0$. Note that the neighboring undershoots around $x=0$ are inherent in the PA operator for $m=3$. Yet, the oscillations left and right from these undershoots stem from the Gibbs phenomenon and thus are parasitical. Applying $\ell^{1}$ regularization to enhance sparsity of the PA transform, however, is able to correct these spurious oscillations, resulting in a sparse representation of the PA transformation in Figure 3(d). The corresponding nodal values, illustrated in Figure 3(c), now approximate the nodal values of the true solution accurately.

For the numerical test in Figure 3, we have chosen the parameters $K=400$, $\mu=0.005$ ($\lambda=4\cdot 10^{2}$), $\beta=20$, $\alpha=0.0001$, and $tol=0.001$. The Lagrangian multipliers have been initialized with $\underline{\sigma}_{0}=\underline{0}$ and $\underline{\delta}_{0}=\underline{0}$.

3.6 Preservation of mass conservation

An essential property of standard (DG) finite element methods for hyperbolic conservation laws is that they are conservative, i.e., that

[TABLE]

holds [69, 3, 35]. Any reasonable shock capturing procedure for hyperbolic conservation laws should preserve this property as well. In particular, in a troubled element $\Omega_{i}$ in which $\ell^{1}$ regularization is applied,

[TABLE]

should hold in order for $\ell^{1}$ regularization to preserve mass conservation of the underlying (DG) method. Unfortunately, (54) is violated when $\ell^{1}$ regularization is applied too naively. This is demonstrated in Figure 4 for a simple discontinuous example $u(x)=\operatorname{sign}(x-0.5)+1$ on $\Omega_{\mathrm{ref}}=[-1,1]$ with mass $\int_{\Omega_{\mathrm{ref}}}u\mathrm{d}x=1$.

The (red) crosses illustrate the absolute difference between the mass of the polynomial interpolation $u_{p}$ and its sparse reconstruction $u_{p}^{\mathrm{spar}}$, i.e.

[TABLE]

for increasing polynomial degrees $p=4,\dots,20$. We observe that a naive application of $\ell^{1}$ regularization might destroy mass conservation. Thus, in the following, we present a simple fix for this problem.

Remark 4**.**

A generic approach — also to, for instance, ensure TVD and entropy conditions — is to add additional constraints to the optimization problem (26), e.g.,

[TABLE]

for conservation of mass to be preserved. In the present case, where $u_{p},v\in\mathbb{P}_{p}\left(\Omega_{\mathrm{ref}}\right)$ are expressed with respect to basis functions $\{\varphi_{k}\}_{k=0}^{p}$ with zero average for $k>0$, the conditions in (56) are easily met. Specifically,

[TABLE]

with

[TABLE]

implies that the additional constraint in (56) reduces to

[TABLE]

Basis functions $\{\varphi_{k}\}_{k=0}^{p}$ with zero average for ${k>0}$ are, for instance, given by the orthogonal basis (OGB) of Legendre polynomials.

We now propose the following simple algorithm to repair mass conservation in the $\ell^{1}$ regularization:

The advantage of this additional step is demonstrated in Figure 4 as well, where the absolute difference between the mass of $u_{p}$ and of its sparse reconstruction with additional mass correction, denoted by $u_{p}^{\mathrm{spar,corr}}$, is illustrated by (black) squares. In contrast to the sparse reconstruction without mass correction, illustrated by (red) crosses, $u_{p}^{\mathrm{spar,corr}}$ is demonstrated to preserve mass nearly up to machine precision ($\approx 10^{-16}$). Finally, we note that for the test illustrated in Figure 3(d), $\ell^{1}$ regularization with and without mass correction resulted in the same approximations, due to $u$ being an odd function. Thus, we omit those results. We present a flowchart in Figure 5 illustrating the proposed procedure for a fixed time $t<t_{\mathrm{end}}$.

4 Numerical tests

We now numerically demonstrate the application of $\ell^{1}$ regularization (with and without mass correction) to a nodal DG method for the inviscid Burgers’ equation, the linear advection equation, and a nonlinear system of conservation laws. Our results show that $\ell^{1}$ regularization provides increased accuracy of the numerical solutions. In all numerical tests, we use a PA operator of third order and choose the same parameters as before, i.e., $\lambda_{\text{max}}=4\cdot 10^{2}$, $K=400$, $\beta=20$, $\alpha=0.0001$, and $tol=0.001$. We have made no effort to optimize these parameters.

4.1 Inviscid Burgers’ equation

We start our numerical investigation by considering the inviscid Burgers’ equation

[TABLE]

on $\Omega=[0,2]$ with initial condition

[TABLE]

and periodic boundary conditions. For this test case, a shock discontinuity develops in the solution at $x=1$.

In the subsequent numerical tests, the usual local Lax–Friedrichs flux

[TABLE]

with maximum characteristic speed $\alpha_{\mathrm{max}}=\max\{|u_{+}|,|u_{-}|\}$ is utilized. Further, we use the third order explicit strong stability preserving Runge–Kutta method with three stages (SSPRK(3,3)) for time integration (see [70]) and choose the ramp parameter ${\kappa=0.8}$ in (31).

Figure 6 illustrates the numerical solutions of the above test problem at time $t=0.345$. In Figure 6(a) a polynomial degree of $p=4$ was used on $I=7$ equidistant elements, while in Figure 6(b) a polynomial degree of $p=10$ was used on $I=7$ equidistant elements and in Figure 6(c) a polynomial degree of $p=18$ was used on $I=7$ equidistant elements. All tests show spurious oscillations for the numerical solution without $\ell^{1}$ regularization. These oscillations are significantly reduced when $\ell^{1}$ regularization is applied to the numerical solution after every time step. As a consequence, Figures 6(d), 6(e), 6(f) illustrate how the pointwise error of the numerical solutions utilizing $\ell^{1}$ regularization (with and without mass correction) is reduced compared to the numerical solution without $\ell^{1}$ regularization. Again, considering an odd function, only very slight differences are observed between $\ell^{1}$ regularization with and without mass correction. The reference solution $u$ was computed using characteristic tracing.

We now extend the above error analysis for the $\ell^{1}$ regularization. Table 2 lists the different common types of errors of the numerical solution for the above test problem (60), (61) and a varying polynomial degree $p$ as well as a varying number of equidistant elements $I$. Here we consider the global errors with respect to the discrete $\underline{\underline{M}}\,$-norm, as defined in (8), approximating the continuous $L^{2}$-norm and the discrete $1$-norm, given by $\mathinner{\!\left\lVert u_{p}\right\rVert}_{1}=\sum_{k=0}^{p}\omega_{k}\left|u_{p}(x_{k})\right|=\left\langle{\underline{1},\,\left|\underline{u}\right|}\right\rangle_{\underline{\underline{M}}\,}$ on the reference element. The discrete $\infty$-norm is given by $\mathinner{\!\left\lVert u_{p}\right\rVert}_{\infty}=\max_{k=0,\dots,p}|u_{p}(x_{k})|$. For the discrete $\underline{\underline{M}}\,$-norm and discrete $1$-norm, the global norms are defined by summing up over the weighted local norms, i.e.,

[TABLE]

where $u_{p}^{(i)}$ denotes the numerical solution (polynomial approximation) on the $i$th element $\Omega_{i}$.

Table 2 demonstrates that for almost all these norms as well as combinations of polynomial degrees $p=4,5,6,7,8,9$ and number of elements $I=15,32,63,127$, the numerical solution with $\ell^{1}$ regularization is more accurate than the numerical solution without $\ell^{1}$ regularization. Further, we observe just a slight difference in accuracy for $\ell^{1}$ regularization with and without mass correction. See, for instance, $p=5$ in Table 2. We only utilize odd number of elements, so that the shock discontinuity arises in the interior of an element. By using an even number of elements, on the other hand, the shock would arise at the interface between two elements and the error analysis of the $\ell^{1}$ regularization would be blurred by the dissipation added by the numerical flux.

In Table 2, all cases where the accuracy is increased or remains the same by applying $\ell^{1}$ regularization are flagged with a checkmark. It should be stressed that no effort was made to optimize the parameters in the $\ell^{1}$ regularization. In particular, the parameters in the $\ell^{1}$ regularization have not been adapted to the specific choice of the polynomial degree $p$, the number of elements $I$, or the number of time steps used in the solver. We think that by further investigations of the parameters in the $\ell^{1}$ regularization, all above errors could be further reduced.

Finally, special attention should be given to the cases flagged by an exclamation mark (!). In these cases, the numerical solver without $\ell^{1}$ regularization broke down completely. Yet, when $\ell^{1}$ regularization was applied, the same computations yielded fairly accurate numerical solutions.

4.2 Linear advection equation

The prior test case featured a shock discontinuity and might not fully reflect the behavior of the $\ell^{1}$ regularized DG method in smooth regions. Note that for the underlying spectral DG method used in this work, given $p+1$ nodes, the optimal order of convergence is $p+1$ in sufficiently smooth regions. Especially in smooth regions, it is desirable that the convergence properties of the underlying method are preserved when using the $\ell^{1}$ regularization method. Convergence of the $\ell^{1}$ regularized DG method, however, also depends on the convergence of the PA operator. In all numerical tests presented in this work, a PA operator of order three, i.e., $L_{3}$, is used. Hence, if falsely activated (false positive; see Remark 3) by the discontinuity sensor (30), $\ell^{1}$ regularization might affect the order of convergence of the underlying method in smooth regions. In our numerical tests, we observed the sensor to be fairly reliable for sufficient resolution, and $\ell^{1}$ regularization did not get activated in smooth regions.

To investigate this potential drawback and demonstrate the reliability of the $\ell^{1}$ regularized DG method in smooth regions, we now consider the linear advection equation

[TABLE]

on $\Omega=[0,2]$ with initial condition

[TABLE]

and periodic boundary conditions, which provides a smooth solution for all times. Table 3 lists comparative errors for DG with and without the $\ell^{1}$ regularization and with the additional mass conservation correction term at time $T=2$.

Our results indicate that $\ell^{1}$ regularization is only activated, and thus affects convergence of the underlying method, if the solution is heavily underresolved. This can be noted in Table 3 from $p=3$ and $N=2,4$ as well as $p=4,5,6$ and $N=2$. For $p=7$, even $N=2$ provides sufficient resolution and the $\ell^{1}$ regularization is not activated. For all numerical solutions, using at least $N=8$ elements, $\ell^{1}$ regularization does not affect accuracy and convergence of the underlying method in smooth regions. Finally, we note that in the cases where the numerical solution is heavily underresolved and $\ell^{1}$ regularization is activated, $\ell^{1}$ regularization with additional mass correction (l1-mc) provides slightly more accurate solutions than $\ell^{1}$ regularization without mass correction (l1). Finally, we note that because the $\ell^{1}$ regularization is only activated in elements containing discontinuities, the efficiency of our new method is comparable to the underlying DG method. Moreover, as in [24], we observed that for our specific set of test problems we were able to maintain stability for time step sizes larger than the standard CFL constraints suggest. Theoretical justification for this will be part of future investigations.

4.3 Systems of conservation laws

We now extend our hybrid $\ell^{1}$ regularized DG method to the nonlinear system of conservation laws

[TABLE]

in the domain $\Omega=[0,2]$. System (66) originates from a truncated polynomial chaos approach for Burgers’ equation with uncertain initial condition [71, 72, 73]. In this context, $u_{0}$ models the expected value of the numerical solution while $u_{1}^{2}$ approximates the variance. For the spatial semidiscretization of (66) we follow [73, 10], where a skew-symmetric formulation

[TABLE]

was proposed. Here $\langle\varphi_{i}\varphi_{j}\varphi_{k}\rangle$ denotes the triple product $\int\varphi_{i}(\xi)\varphi_{j}(\xi)\varphi_{k}(\xi)\omega(\xi)\operatorname{d\!}\xi$, $\circ$ denotes the componentwise (Hadamard) product of two vectors, and $\{\varphi_{k}\}$ is a set of orthogonal polynomials (typically Hermite polynomials are used).

It was further proved in [73] that (67) yields an entropy conservative semidiscretization when combined with the entropy conservative flux $\underline{f_{k}}^{\mathrm{num}}$ presented in [73]. An entropy stable semidiscretization is thus obtained by adding a dissipative term $-\underline{\underline{Q}}\,(\underline{u_{k}}_{+}-\underline{u_{k}}_{-})$ to the entropy conservative flux. Here we simply use a local Lax–Friedrichs type dissipation matrix

[TABLE]

where $\mathinner{\!\left\lvert{\alpha}(\pm)\right\rvert}$ is the largest absolute value of all eigenvalues of the Jacobian $f^{\prime}(u(\pm))$.

Even though the skew-symmetric formulation (67) combined with an appropriate numerical flux yields an entropy stable scheme, this test case demonstrates that additional regularization is still necessary to obtain reasonable numerical solutions. This was, for instance, stressed in [10]. Here we demonstrate how $\ell^{1}$ regularization enhances the numerical solution for the nonlinear system (66) with periodic boundary conditions and initial condition

[TABLE]

where $r=0.5$.

For the more general case of systems of conservation laws, we propose a straightforward extension of our $\ell^{1}$ regularization technique. Specifically, the PA sensor (30) is applied to every conserved variable $u_{k}$ separately and, once a discontinuity is detected, $\ell^{1}$ regularization (27) is performed for the respective variable. For this test case, the ramp parameter in (31) has been chosen as $\kappa=0.9$.

Figure 7 illustrates the results of $\ell^{1}$ regularization (with and without mass correction) for the above described test case and for an entropy stable numerical flux. In all subsequent tests $I=100$ equidistant elements and a polynomial basis of degree $p=6$ have been used. Further, for the $\ell^{1}$ regularization, the same parameters as before have been used, i.e., $\lambda_{\text{max}}=4\cdot 10^{2}$, $K=400$, $\beta=20$, $\alpha=0.0001$, and $tol=0.001$.

Note that while the numerical solution without $\ell^{1}$ regularization shows heavy oscillations in both components, the numerical solution with $\ell^{1}$ regularization provides a significantly sharper profile. Further, by consulting Figures 7(d) and 7(e), it should be stressed that only $\ell^{1}$ regularization with additional mass correction is able to capture the exact shock location. Due to missing conservation, $\ell^{1}$ regularization without mass correction results in a slightly wrong location for the shock. Finally, Figures 7(c) and 7(f) illustrate the energy of the different methods over time. We note from these figures that $\ell^{1}$ regularization (with and without mass correction) slightly increase the energy in this test case.

In order to further emphasize the effect of $\ell^{1}$ regularization (with and without mass correction), similar results using an entropy conservative numerical flux at the interfaces between elements are shown in Figure 8. Once again, $\ell^{1}$ regularization is demonstrated to improve the numerical solution. The best results are obtained when $\ell^{1}$ regularization is combined with mass correction. In particular, only $\ell^{1}$ regularization with additional mass correction is able to accurately capture the shock at the right location. Finally, consulting Figures 8(c) and 8(f), $\ell^{1}$ regularization decreases the energy overall. Yet, it is demonstrated once more that an entropy inequality is not satisfied by $\ell^{1}$ regularization, i.e., the energy might increase as well as decrease by utilizing $\ell^{1}$ regularization. Thus, future work will focus on incorporating energy stability (as well as other properties like TVD or positivity) by additional constraints in the minimization problem (26).

5 Concluding remarks

We have presented a novel approach to shock capturing by $\ell^{1}$ regularization using SE approximations. Our work not only is distinguished from previous studies [20, 21, 22, 23, 25] by focusing on discontinuous solutions but further by promoting sparsity of the jump function instead of the numerical solution itself. By approximating the jump function with the (high order) PA operator, we help to eliminate the staircase effect that arises for classical TV operators. Our results demonstrate that it is possible to efficiently implement a method that yields increased accuracy and better resolves (shock) discontinuities. In particular, no additional time step restrictions are introduced, in contrast to artificial viscosity methods when no care is taken in their construction. This approach for solving numerical conservation laws was first used in [24], where the Lax–Wendroff scheme and Chebyshev and Fourier spectral methods were used as the numerical PDE solver. Our method improves upon the approach in [24] in two ways. First, we employ the SE approximation for solving the conservation law, which allows element-to-element variations in the optimization problem. In particular, $\ell^{1}$ regularization is only activated in troubled elements, which enhances accuracy and efficiency of the method. Second, in the process we proposed a novel discontinuity sensor based on PA operators of increasing orders, which is able to flag troubled elements as well as to steer the amount of regularization introduced by the sparse reconstruction.

Numerical tests demonstrate the method using a nodal collocation-type discontinuous Galerkin method for the inviscid Burgers’ equation, the linear advection equation, and a nonlinear system of conservation laws. Our results show that the method yields improved accuracy and robustness.

No effort was made in our study to optimize any of the parameters involved in solving the optimization problem. This will be addressed in future work, along with the possibility to include additional constraints (e.g., for entropy, TVD, and positivity constraints), since preliminary results presented here are encouraging. The generalization of the approach itself to higher dimensions is straightforward and has already been demonstrated in [24]. Of interest, however, would be the extension of the proposed approach to other classes of methods, such as finite volume methods. We believe $\ell^{1}$ regularization might be an important ingredient to make high order methods viable in several research applications.

Acknowledgements

The authors would like to thank Chi-Wang Shu (Brown University) for helpful advice. Further, the authors would like to thank the anonymous referees for many helpful suggestions, resulting in an improved presentation of this work. Jan Glaubitz’ work was supported by the German Research Foundation (DFG, Deutsche Forschungsgemeinschaft) under grant SO 363/15-1. Anne Gelb’s work was partially supported by AFOSR9550-18-1-0316 and NSF-DMS 1502640.