Shock Regularization with Smoothness-Increasing Accuracy-Conserving Dirac-Delta Polynomial Kernels
Wissink B.W, Jacobs G.B., Ryan J.K., Don W-S, van der Weide E.T.A
SShock Regularization with Smoothness-IncreasingAccuracy-Conserving Dirac-Delta Polynomial Kernels
B.W. Wissink a,d , G.B. Jacobs a , J.K. Ryan b ,W.S. Don c and E.T.A. van der Weide da Department of Aerospace Engineering, San Diego StateUniversity, San Diego, CA b School of Mathematics, University of East Anglia,Norwich, UK c School of Mathematical Sciences, Ocean Universityof China, Qingdao, China d Faculty of Engineering Technology, University of Twente,Enschede, The Netherlands corresponding author, [email protected] 25, 2017
Abstract
A smoothness-increasing accuracy conserving filtering approach to the regularization ofdiscontinuities is presented for single domain spectral collocation approximations of hyperbolicconservation laws. The filter is based on convolution of a polynomial kernel that approximatesa delta-sequence. The kernel combines a k th order smoothness with an arbitrary numberof m zero moments. The zero moments ensure a m th order accurate approximation of thedelta-sequence to the delta function. Through exact quadrature the projection error of thepolynomial kernel on the spectral basis is ensured to be less than the moment error. A numberof test cases on the advection equation, Burger’s equation and Euler equations in 1D and 2Dshown that the filter regularizes discontinuities while preserving high-order resolution. Keywords
Shock capturing, Hyperbolic conservation laws, Regularization, Dirac-Delta, Chebyshev collocation,Filtering 1 a r X i v : . [ m a t h . NA ] A ug Introduction
Shock capturing with high-order spectral methods is well known to be plagued by Gibbs phenom-ena in the solution. Nonlinear polynomial reconstruction schemes such as the nonlinear weightedessentially non-oscillatory (WENO) finite difference schemes on a uniformly spaced grid that havebeen very successful [1] do not extend well to global polynomial based Chebyshev and Legendrecollocation (spectral) methods. Since spectral methods rely on high-order global basis functions,slope limiters are often applied to suppress Gibb’s oscillations. The most common approach is touse explicit Runge Kutta Discontinuous Galerkin (RKDG) methods with min-mod slope limitersintroduced by Cockburn [2], [3]. A cost-effective alternative to limiting is the artificial viscosity(AV) approach, e.g. [4], [5], [6] where artificial higher even order differential terms are added tothe equations to dissipate the high frequency waves or smoothen the small scale structures. Whilethis approach is very stable and more accurate than lower order alternatives, there is no formalproof of higher-order resolution and accuracy. Yet another approach is to use filtering, which hassuccessfully employed in simulating shocked flow [7, 4]. The high order filter used there is thevariable order exponential filter, which does not satisfy all the criteria for the definition of filter aslaid out in [8] exactly but asymptotically.A yet to be explored technique for shock capturing is the use of SIAC-like filters. The typical appli-cation of SIAC filtering is to obtain superconvergence. This is accomplished by using informationthat is already contained in the numerical solution to increase the smoothness of the field and toreduce the magnitude of the errors. The solution, as a post-processing step, is convolved against aspecifically designed kernel function once at the final time. The foundations for this postprocessorwere established by Bramble and Schatz [9]. They showed that the accuracy of Ritz-Galerkin dis-cretizations can be doubled by convolving the solution against a certain kernel function. Cockburnet al. [10] used the ideas of Bramble and Schatz and those of Mock and Lax [11] to demonstrate thatthe postprocessor is also suitable for DG schemes. They proved that, for a certain class of linearhyperbolic equations with sufficiently smooth solutions, the postprocessor enhances the accuracyfrom order k + 1 to order 2 k + 1 in the L -norm, where k is the polynomial degree of the originalDG approximation. This postprocessor relies on a symmetric convolution kernel consisting of 2 k + 1B-splines of order k + 1.In [12] a regularization technique is developed for the Dirac- δ source terms in hyperbolic equationsthat is an excellent candidate for a kernel of a SIAC-like regularization filter of shock discontinu-ities. The technique is based on a class of high-order compactly supported piecewise polynomialsintroduced in [13]. The piecewise polynomials provide a high-order approximation to the Dirac-Delta whose overall accuracy is controlled by two conditions: the number of vanishing momentsand smoothness. SIAC kernels have similar smoothness and moment properties as the Dirac-Delta kernel, but are based on piecewise continuous B-splines instead of polynomials. In [14] it wasshown that SIAC-like filters based on the compactly supported Dirac-Delta kernels capture particle-fluid interface discontinuities with higher-order resolution in single domain spectral solutions. Thecompactness of support is closely related to the moment condition. Smoothness properties yieldhigher-order resolution away from the discontinuity.2n the present work we have developed SIAC-like filters based on the high-order Dirac-Delta kernelfor the regularization of shocks and discontinuities. The filter operation is based on convolution ofthe solution with the high-order Dirac-Delta kernels. The formulation of the operation is derived andwritten in matrix-vector multiplication form to allow for an efficient and simple implementation. Thefilters are tested on a spectral solver for hyperbolic conservation laws, such as the one dimensionalscalar linear advection equation and scalar nonlinear Burgers’ equation, and both one- and two-dimensional nonlinear Euler equations. We show that, for the case of the linear advection equation,filtering discontinuous initial conditions yields a stable, high-order resolution and with the designedorder of accuracy away from the discontinuity.In the area where the solution is the advected (filtered) initial condition, the designed order ofaccuracy is achieved and depends on the support width and number of vanishing moments on thekernel. The support width is made proportional to grid-spacing and also depends on the numberof vanishing moments.In the case of the non-linear Burgers’ and Euler equations, the filter has to be applied at every timestep in order to maintain stability by smoothing shocks and discontinuities during the simulation.This leads to a summation of filter-errors and smearing of the solution near discontinuities. Weshow that by choosing a sufficiently accurate Dirac-Delta kernel and small support-width, theseeffects are minimized and design order of accuracy is obtained away from the discontinuity as well.In Section 2 the filter-operation is derived and background information about SIAC filters, thehigh-order Dirac-Delta functions and the Chebyshev collocation method is provided. Next, Section3 presents and discusses the numerical results. Finally, Section 4 summarizes the results and givesan outlook for future work. Polynomial based spectral methods, in this particular case, the Chebyshev collocation method, arecommonly used in the discretization of spatial derivatives in PDE’s since the order of convergence,for a sufficiently smooth function, depends only on the smoothness of the function, also known asspectral accuracy. For example, the spectral approximation of a C p function is at least O ( N − p ).In the case of an analytical function, the order of convergence of the approximation is exponential[15]. In the following the Chebyshev collocation method is briefly described for the purpose ofintroducing notation. For an overview on spectral methods we refer to [15].The collocation method is based on polynomial interpolation of a function u ( x ), and can be ex-pressed as u N ( x ) = N (cid:88) j =0 u ( x j ) l j ( x ) , l j ( x ) = N (cid:89) k =0 ,k (cid:54) = j x − x k x j − x k , j = 0 , . . . , N, (1)3here x j , j = 0 , . . . , N are the collocation points and l j ( x ) are the Lagrange interpolation polyno-mials of degree N . To determine the derivative of the function u ( x ) at the collocation points x i , u (cid:48) ( x i ), one can simply take the derivative of the Lagrange interpolating polynomial as ∂u ( x i ) ∂x ≈ N (cid:88) j =0 u ( x j ) l (cid:48) j ( x i ) . (2)or, written compactly in the matrix-vector multiplication form as −→ u (cid:48) = D −→ u , (3)where the differentiation matrix D i,j = l (cid:48) j ( x i ). In the case of the Chebyshev collocation method,the collocation points are x i = − cos( iπ/N ) , i = 0 , . . . , N. (4)To integrate the resulting system of ordinary differential equations (ODE) in time, we employ thethird order Total Variation Diminishing (TVD) Runge-Kutta scheme u (1) = u n + ∆ tL ( u n ) u (2) = 34 u n + 14 u (1) + 14 ∆ tL ( u (1) ) (5) u n +1 = 13 u n + 23 u (2) + 23 ∆ tL ( u (2) ) . In [12], a regularization technique based on a class of high-order, compactly supported piecewisepolynomials is developed that regularizes the time-dependent, singular Dirac-Delta sources in spec-tral approximations of hyperbolic conservation laws. This regularization technique provides higher-order accuracy away from the singularity.For the purpose of clarity in the following discussion, we shall define the compact support domainas Ω ε = [ ε − x, x + ε ] and Ω εi = [ x i − ε, x i + ε ] centered at x = x i with a given support width ε > δ m,kε ( x ), that is an approximationof the Dirac-Delta function δ ( x ), δ m,kε ( x ) = (cid:40) ε P m,k (cid:0) xε (cid:1) x ∈ Ω ε , , (6)where ε > m , and the number of continuousderivatives at the endpoints of the compact support k , respectively. The M = m + 2( k + 1) degree4olynomial P m,k ( x ) is uniquely determined by the following properties (cid:0) P m,k (cid:1) ( i ) ( ±
1) = 0 , i = 1 , ..., k, (7) (cid:90) − P m,k ( ξ ) dξ = 1 , (8) (cid:90) − ξ i P m,k ( ξ ) dξ = 0 , i = 1 , ..., m, (9)in which (7) determines the number of continuous derivatives at the endpoints ( k ), (8) states that thefunction integrates to unity as a Dirac-Delta function, and (9) determines the number of vanishingmoments ( m ). In Fig. 1, we show the polynomial approximation of the Dirac-Delta kernel δ m,kε ( x )with ( m, k ) = (3 ,
8) and ( m, k ) = (5 ,
8) with scaling parameter ε = 1. It has been shown that theDirac-Delta approximation δ m,kε ( x ) has an accuracy of O ( ε m +1 ).The procedure for the generation of the polynomials P m,k ( ξ ) is described in [13]. A few examplesfor low values of m and k are given below. • The polynomials with one vanishing moment m = 1 and with k = 0 , P , = 34 (1 − ξ ) , P , = 1516 (1 − ξ + ξ ) , P , = 3532 (1 − ξ + 3 ξ − ξ ) . (10) • The polynomials with three vanishing moments m = 3 and with k = 0 , P , = 1532 (3 − ξ + 7 ξ ) , P , = 10564 (1 − ξ + 7 ξ − ξ ) , (11) P , = 315512 (3 − ξ + 42 ξ − ξ + 11 ξ ) . (12) The Smoothness-Increasing Accuracy-Conserving (SIAC) filter was developed in the context ofsuperconvergence extraction and error reduction. SIAC filtering [16, 17, 18, 19] exploits the idea ofsuperconvergence in the underlying method to reduce oscillations in the errors, reduce errors andincrease convergence rates. SIAC has its basis in the work by Bramble and Schatz [9] and Cockburn,Luskin, Shu and S¨uli [10]. It has been extended to include boundary filtering and nonuniform meshes[20, 21, 22, 23, 24] as well as to increase computational efficiency [25, 26]. It has proven effectivefor applications in visualization [27, 28] and it is expected that this superconvergence extractiontechnique can act as a natural filter for EFLES.The traditional application of SIAC filtering takes the numerical approximation, u h ( x ), and con-volves it against a specially designed kernel, u (cid:63)h ( x ) = 1 H (cid:90) ∞−∞ K m +1 ,(cid:96) (cid:18) y − xH (cid:19) u h ( y ) dy, (13)5igure 1: Dirac-Delta kernel δ m,kε for ( m, k ) = (3 ,
8) and ( m, k ) = (5 ,
8) with scaling parameter ε = 1.where K m +1 ,(cid:96) is a linear combination of m B-splines of order (cid:96) and H is a scaling typically relatedto a mesh quantity.The symmetric form of the post-processing kernel can be written as K m +1 ,(cid:96) ( x ) = (cid:88) γ c m +1 ,(cid:96)γ ψ ( (cid:96) ) ( x − γ ) , (14)where ψ ( (cid:96) ) ( x ) is obtained by convolving the characteristic function over the interval ( − , ) withitself (cid:96) + 1 times and the coefficients c m +1 ,(cid:96)γ ∈ R . The form of the SIAC kernel is similar to the form of the mixed polynomial [12, 13]. Primaryproperties of the SIAC that make the SIAC kernel suitable for regularization include: • ψ ( (cid:96) ) ( x ) can be expressed as a linear combination of Delta-functions, using the property d α ψ ( (cid:96) ) ( x ) dx α = ∂ αH ψ ( (cid:96) − α ) ( x ) , where ∂ H represents a divided difference. • The SIAC kernel can reproduce polynomials of degree m. This is equivalent to the conditions – (cid:82) R K m +1 ,(cid:96) ( ξ ) dξ = 1 . – (cid:82) R ξ i K m +1 ,(cid:96) ( ξ ) dξ = 0 for i = 1 , . . . , m. These are similar to the conditions on mixed polynomials. However, unlike the mixed polyno-mials, the SIAC kernel does not require the smoothness of the kernel even though it vanishes,at the endpoints of the compact support.
In [14], an extension from [12], singular source terms are expressed as weighted summations ofDirac-Delta kernels that are regularized through approximation theory with convolution operators.The regularization is obtained by convolution with the high-order compactly supported Dirac-Delta6ernel, (6), whose overall accuracy is controlled by the number of vanishing moments m , degreeof smoothness k and length of the support ε . In this work, the Dirac-Delta kernel is used toregularize time dependent discontinuous solutions . Suppose the solution is given by the variable u ( x ) , x ∈ [ − , u ( x ), follows from the convolution of u ( x ) with theDirac-Delta kernel as ˜ u ( x ) = (cid:90) − u ( τ ) δ m,kε ( x − τ ) dτ. (15)or, simply, ˜ u ( x ) = (cid:90) Ω ε u ( τ ) δ m,kε ( x − τ ) dτ, (16)since the Dirac-Delta kernel is zero outside its compact support Ω ε .To apply the filtering operation, we need to choose the number of vanishing moments m , the numberof continuous derivatives k and the scaling parameter ε .We make the following notes: • The number of vanishing moments and support width determine the accuracy of the Delta-kernel and therefore the error introduced by the filtering operation in smooth regions ofthe solution. In [14] it was proven that the filtering error is O ( ε m +1 ), provided the scalingparameter is ε = O ( N ( − k/m + k +2) ). The requirement on the scaling parameter follows fromthe fact that the error in the quadrature rule used to evaluate the convolution integral hasto be smaller than the O ( ε m +1 ) accuracy of the Dirac-Delta approximation. More vanishingmoments reduces the filtering error, however it requires a wider support, leading to a widerregularization zone. • The smoothness of the Dirac-Delta kernel is controlled by the number of continuous derivatives k at the endpoints. In [14] and [12] it is shown that k controls the smoothness of the transitionbetween the regularized source and the solution and thus controls the order of convergenceaway from the source. When filtering the entire solution there is no such transition and thusthe influence of k is minor. This is confirmed by numerical experiments.Extension to two dimensions follows straightforwardly from the tensor product of the one-dimensionalDelta-function: δ m,kε ( x, y ) = δ m,kε ( x ) ⊗ δ m,kε ( y ) (17)Fig. 2 shows the two-dimensional equivalents of the one-dimensional kernels shown in Fig. 1.As in the one-dimensional case, the filtering is based on the convolution of the solution with theDelta-kernel, ˜ u ( x, y ) = (cid:90) Ω εx (cid:90) Ω εy u ( τ, η ) δ m,kε ( x − τ, y − η ) dτ dη. (18)7 m, k ) = (3 ,
8) ( m, k ) = (5 , m, k ) = (3 ,
8) and (Right) ( m, k ) = (5 ,
8) with scaling parameter ε = 1. Filtering of the interpolant of the solution (1), leads to,˜ u N ( x ) = (cid:90) Ω ε (cid:34) N (cid:88) i =0 u ( x i ) l i ( τ ) (cid:35) δ m,kε ( x − τ ) dτ = N (cid:88) i =0 u ( x i ) S i ( x ) , (19)after interchanging the summation and integration, where the filtering function S i ( x ) is S i ( x ) = (cid:90) Ω ε l i ( τ ) δ m,kε ( x − τ ) dτ. (20)Hence, one has, at the collocation points, −→ ˜ u = S −→ u , (21)where the ( N + 1) × ( N + 1) filtering matrix S has elements S i,j = (cid:90) Ω εj l i ( τ ) δ m,kε ( x j − τ ) dτ. (22)In two dimensions, the filtering operation can be written compactly as˜ U = S x US T y , (23)where S x and S y are the one-dimensional filtering matrix in x − and y − direction respectively andthe superscript T denotes transpose. Remark 1
Near the boundaries, the compact support of the high-order Dirac-Delta function, δ m,kε ,extends out of the domain and hence, the data can not be filtered in that case. This means that ˜ u ( x ) = u ( x ) for | x | > − ε . The filtering matrix S can be precomputed and stored for later use aslong as the filter parameters remain unchanged. .6 Clenshaw-Curtis quadrature The one remaining important issue that needs to be addressed is how to evaluate the integrals inEq. 19. Since the high-order Dirac-Delta function is a polynomial of degree M = m + 2( k + 1) andthe Lagrange interpolation polynomial is of degree N −
1, the integrand is a polynomial of degree M + N − N an appropriate quadrature rule is preferred. For this purpose Clenshaw-Curtis quadrature will beused, that is, (cid:90) Ω εi l n ( τ ) δ m,kε ( x i − τ ) dτ ≈ Q (cid:88) q =0 w q l n ( x q ) δ m,kε ( x i − x q ) . (24)where Q is the number of Chebyshev Gauss-Lobatto quadrature nodes used in the compact supportdomain Ω εi , that is, x q = x i − ε cos (cid:18) πqQ (cid:19) , q = 0 , . . . , Q, (25)and w q are the corresponding weights. Clenshaw Curtis quadrature exactly evaluates polynomialsof degree Q −
1. Hence, if one takes Q = M + N , the integrals are evaluated exactly. Also, theweights w q can be precomputed using fast Fourier transform (FFT). The theoretical requirement for the scaling parameter, ε = O ( N ( − k/m + k +2) ), to assure O ( ε m +1 )accuracy in the convolution operation is based on the requirement that the error in the quadraturerule has to be smaller than the error in the Dirac-Delta approximation [14]. Since the convolutionintegrals in this work are solved exactly using Clenshaw-Curtis quadrature, this criterium can berelaxed. For shock capturing we choose the scaling parameter to be proportional to the grid spacingand to guarantee that a minimal of at least two neighboring collocation points will be located insidethe compact support width Ω εi for any point x = x i . In this case, the scaling parameter will beexpressed in terms of the number of points N d the kernel spans at the center of the domain x = x N/ ,assuming N is even, that is, ε = sin ( πN d / (2 N )) . (26)The value of N d is determined empirically and chosen to ensure stable converging results. It dependson the number of vanishing moments m of the Dirac-Delta kernel and whether the filter is appliedin a linear or non-linear equation. For implementation in the linear advection equation, only theinitial condition is filtered while for a nonlinear PDE’s, such as the Burgers’ equation and the Eulerequations, the solution is filtered at every Runge Kutta time step. This is because of the formationof a finite space-time singularity by the nonlinearity of the equations. However, filtering can leadto smearing of the discontinuity and the summation of filter errors in smooth regions. In order tolimit these less desirable effects, we use a smaller scaling parameter, ε , as compared to the linearadvection equation. 9 emark 2 Near the boundary the filter cannot be convoluted due to the symmetry of the convolutionkernel. In this paper those few points are reset to be the exact solution for the nonlinear PDEs inorder to isolate and study solely the effects of the Dirac-Delta kernels. Furthermore, to avoid theeffect of variable time step ∆ t , we fixed the stable time step of ∆ t = 1 × − in all the simulationsperformed below in order to be able to compare results. We shall consider the simple one-dimensional linear advection equation with a discontinuity asinitial condition on the domain − < x < ∂u∂t + ∂u∂x = 0 , (27) u ( x,
0) = (cid:26) sin( πx ) − . x ≤ − . πx ) + 0 . x > − . ,u ( − , t ) = sin( π ( − − t )) − . . The initial condition is first filtered by the filtering matrix S . For the results below we used aDirac-Delta kernel with ( m, k ) = (3 ,
8) and N d = 13.Figure 3: Linear advection equation: Analytical and spectral solution (Left), and pointwise error(Right) at time t = 1, for the five different grids with ( m, k ) = (3 ,
8) and N d = 13.In Fig. 3, the point-wise error clearly shows the advected filter error for x > x = 0 .
28 and follows the theoretical value of O ( ε m +1 ).10igure 4: Linear advection equation: Detail view of regularization zone (Left) and error convergenceat x = 0 .
28 compared to the theoretical value O ( ε m +1 ) (Right) for the five different grids with( m, k ) = (3 ,
8) and N d = 13. Next, we consider the Burgers’ equation, ∂u∂t + 12 ∂u ∂x = 0 , (28) u ( x,
0) = − sin( πx ) , (29) u ( ± , t ) = 0In this case, the solution is filtered at the end of every Runge-Kutta time step. For the resultsbelow we used a kernel with ( m, k ) = (3 ,
8) and N d = 2 . t = 1, for the five different grids with ( m, k ) = (3 ,
8) and N d = 2 . N d used here,the smearing is is small. The smearing effect increases with increasing support width. We not thatextending the simulation time until t = 5 does not change the error much as compared to error at t = 1. Only the shock is slightly more dissipated. The 2D unsteady Euler Equations in conservative form are given as ∂U∂t + ∂F x ∂x + ∂F y ∂y = 0 , (30)where the conserved quantities U and the flux vectors F x and F y are given by U = ρρuρvρE , F x = ρuρu + pρuv ( ρE + p ) u , F y = ρuρuvρv + p ( ρE + p ) v . (31)The system is closed by assuming a calorically perfect gas, which relates the pressure to the density,velocity and energy. For the two-dimensional case this gives p = ( γ − (cid:18) ρE − ρ ( u + v ) (cid:19) , (32)in which γ is the specific heat ratio, which is 1 . Sod’s shock tube problem is governed by the 1D Euler equations. The initial conditions for Sod’sshock tube problem are ρ L = 1 . ρ R = 0 . ,p L = 1 . p R = 0 . ,u L = 0 . u R = 0 . . where the left state is in the domain with x < x ≥ m, k ) = (3 ,
8) and N d = 2 .
5. The solutions show that both the shock and the contactdiscontinuity are effectively captured using the filter-operation. A low error is observed outside theregularization zone. 12igure 6: Analytical solution, spectral solution (Left) and pointwise error (Right) of the density at t = 0 .
4, for the five different grids for ( m, k ) = (3 ,
8) and N d = 2 . Shu-Osher’s Problem is also governed by the 1D Euler equations with the initial conditions ofthe primitive variables in the left state ( ρ L , P L , u L ) = (27 / , / , (cid:112) / ρ R , P R , u R ) = (1 + 0 . x ) , , m, k ) = (5 ,
8) and N d = 5 . Results, δ , ε , N d = 6 . Figure 7: Analytical solution, spectral solution (Left) and pointwise error (Right) of the density at t = 0 .
18, for the five different grids for ( m, k ) = (5 ,
8) and N d = 6 . t = 0 .
36, for the five different grids for ( m, k ) = (5 ,
8) and N d = 6 . t = 0 . m, k ) = (5 ,
8) and N d = 6 . m = 5vanishing moments ensures high resolution behind the shock. Fig. 9 provides a detail view of theentropy wave (Left) and the regularization zone (Right). As a 2D test case, the explosion problem is considered. This problem is governed by the 2D Eulerequations, (30), (31) and (32). The equations are solved on a 2 . × . x − y plane. The initial condition consists of the region inside of a circle with radius R = 0 . ,
0) and the region outside of the circle. The flow variables are constant in each of these regions14nd are separated by a circular discontinuity at time t = 0. The two constant states are chosen as: ρ in = 1 . ρ out = 0 . , (33) p in = 1 . p out = 0 . , (34) u in = 0 . u out = 0 . , (35) v in = 0 . v out = 0 . , (36)where the subscripts in and out denote values inside and outside the circle, respectively. Theapproach described by Eleuterio [29] is used to serve as an exact solution. For the results below weused a kernel with ( m, k ) = (3 ,
8) and N d = 2 . t = 0 .
25 for ( m, k ) = (3 ,
8) and N d = 2 . t = 0 .
25, for the fivedifferent grids for ( m, k ) = (3 ,
8) and N d = 2 . Conclusions and future work
In this paper, the use of high-order Dirac-Delta function based filters for the regularization of shocksand discontinuities in combination with a global spectral method is investigated. We have shownthat these filters are able to effectively capture shocks while maintaining high resolution in smoothareas of the solution.As a first next step for future work, the filtering operation has to be extended for the use indiscontinuous spectral element methods. Furthermore the development of a way to apply thefiltering near the boundaries is desired. A possible way could be the development of special one-sided polynomial kernels, analog to the one-sided SIAC kernels for boundary filtering [28], [20].In the present work, the developed filters have been applied globally, i.e. the whole domain isaffected by the filter, thus causing (small) errors in smooth areas and smearing of the discontinuity.It would be beneficial if this effect could be limited. This could be accomplished if only the area ofthe discontinuity is filtered. This however requires some means to identify shocks and discontinuitiesin the solution. One possibility is the use of multi-resolution analysis (MR) [30], [31]. MR-analysisis successfully used in hybrid spectral-WENO methods to switch between high-resolution WENOmethods for shocked regions and computationally efficient spectral methods in smooth regions [32].In a follow on paper, we intend to report on the use the filter-error itself to identify shocks anddiscontinuities. The filter-error e = | u − S u | , (37)describes the shape of the solution in the sense that it is small in smooth areas and large in shockedareas. If a function L ( e ) can be derived that is able to determine which areas should be filtered andwhich not, the filtering operation could be improved and implemented according to˜ u = L ( e ) u + (1 − L ( e )) S u . (38)Thus for areas in which L ( e ) = 1 the solution is not filtered whereas if L ( e ) = 0 the solution is fullyfiltered. This would then effectively limit the smearing and the error introduced in smooth areas. References [1] C.W. Shu and S. Osher. Efficient Implementation of Essentially Non-oscillatory Shock-Capturing Schemes.
Journal of Computational Physics , 77:439–471, 1988.[2] B. Cockburn. Devising discontinuous galerkin methods for non-linear hyperbolic conservationlaws.
J. Comp. Apll. Math. , 128:187–204, 2001.[3] B. Cockburn and C.W. Shu. Runge-kutta discontinuous galerkin methods for convection-dominated problems.
Journal of scientific computing , 16:173–261, 2001.[4] A. Chaudhuri, G.B. Jacobs, W.S. Don, H. Abassi, and F. Mashayek. Explicit discontinuousspectral element method with entropy generation based artificial viscosity for shocked viscousflows.
J. Comp. Phys. , 332:99–117, 2017. 165] T.J. Hughes, L. Franca, and M. Mallet. A new finite element formulation for computationalfluid dynamics: I. symmetric forms of the compressible euler and navier-stokes equations andthe second law of thermodynamics.
Computer Methods in Applied Mechanics and Engineering ,54:223–234, 1986.[6] J.L. Guermond, R. Pasquetti, and B. Popov. Entropy viscosity method for nonlinear conser-vation laws.
Journal of Computational Physics , 230:4248–4267, 2011.[7] W.S. Don. Numerical Study of Pseudospectral Methods in Shock Wave Applications.
J.Comput. Phys. , 110:103–111, 1994.[8] H. Vandeven. Family of spectral filters for discontinuous problems.
Journal of Sientific Com-puting , 6(2):159–192, 1991.[9] J.H. Bramble and A.H. Schatz. Higher Order Local Accuracy by Averaging in the FiniteElement Method.
Mathematics of Computation , 31:94–111, 1977.[10] B. Cockburn, M. Luskin, C.W. Shu, and E. Sli. Enhanced Accuracy by Post-processing forFinite Element Methods for Hyperbolic Equations.
Mathematics of Computation , 72:577–606,2003.[11] M.S. Mock and P.D. Lax. The computation of discontinuous solutions of linear hyperbolicequations.
Comm. Pure Appl. Math. , 18:423–430, 1978.[12] J.P. Suarez, G.B. Jacobs, and W.S. Don. A high-order Dirac-delta regularization with optimalscaling in the spectral solution of one-dimensional singular hyperbolic conservation laws.
SIAMJ. Sci. Comp , 36:1831–1849, 2014.[13] A.K. Tornberg. Multi-dimensional quadrature of singular and discontinuous functions.
BIT ,42:644–669, 2002.[14] J.P. Suarez and G.B. Jacobs. Regularization of singularities in the weighted summation ofDirac-delta functions for the spectral solution of hyperbolic conservation laws.
J. Sci. Comp. , preprint at: https://arxiv.org/abs/1611.05510 .[15] Sigal Gottlieb Jan Hesthaven and David Gottlieb. Spectral Methods for Time-DependentProblems. Cambridge University Press, Cambridge, United Kingdom .[16] J.K. Ryan, C.W. Shu, and H.L. Atkins. Extension of a post-processing technique for discontin-uous Galerkin methods for hyperbolic equations with application to an aeroacoustic problem.
SIAM Journal on Scientific Computing , 26:821–843, 2004.[17] J. King, H. Mirzaee, J.K. Ryan, and R.M. Kirby. Smoothness-Increasing Accuracy-Conserving(SIAC) filtering for discontinuous Galerkin solutions: Improved errors versus higher-order ac-curacy.
Journal of Scientific Computing , 53:129–149, 2012.[18] S. Curtis, R.M. Kirby, J.K. Ryan, and C.W. Shu. Post-processing for the discontinuousGalerkin method over non-uniform meshes.
SIAM Journal on Scientific Computing , 30:272–289, 2007. 1719] L. Ji, P. van Slingerland, J.K. Ryan, and C. Vuik. Superconvergent error estimates for aposition-dependent Smoothness-Increasing Accuracy-Conserving filter for DG solutions.
Math-ematics of Computation , 83:2239–2262, 2014.[20] J.K. Ryan and C.W. Shu. On a one-sided post-processing technique for the discontinuousGalerkin methods.
Methods Appl. Anal. , 10:295–307, 2003.[21] P. van Slingerland, J.K. Ryan, and C. Vuik. Position-dependent Smoothness-IncreasingAccuracy-Conserving (SIAC) filtering for accuracy for improving discontinuous Galerkin solu-tions.
SIAM J. Sci. Comp. , 33:802–825, 2011.[22] X. Li, R.M. Kirby, and C. Vuik. Computationally efficient position-dependent Smoothness-Increasing Accuracy-Conserving (SIAC) filtering: The uniform mesh case.
Preprint , 2013.[23] X. Li, R.M. Kirby, and C. Vuik. Computationally efficient position-dependent Smoothness-Increasing Accuracy-Conserving (SIAC) filtering: The nonuniform mesh case.
Preprint , 2014.[24] H. Mirzaee, L. Ji, J.K. Ryan, and R.M. Kirby. Smoothness-Increasing Accuracy-Conserving(SIAC) post-processing for discontinuous Galerkin solutions over structured triangular meshes.
SIAM Journal on Numerical Analysis , 49:1899–1920, 2011.[25] H. Mirzaee, J.K. Ryan, and R.M. Kirby. Quantification of errors introduced in the numeri-cal approximation and implementation of Smoothness-Increasing Accuracy-Conserving (SIAC)filtering of discontinuous Galerkin (DG) fields.
Journal of Scientific Computing , 45:447–470,2010.[26] H. Mirzaee, J.K. Ryan, and R.M. Kirby. Efficient implementation of Smoothness-IncreasingAccuracy-Conserving (SIAC) filters for discontinuous Galerkin solutions.
Journal of ScientificComputing , 52:85–112, 2012.[27] M. Steffan, S. Curtis, R.M. Kirby, and J.K. Ryan. Investigation of Smoothness-IncreasingAccuracy-Conserving Filters for Improving Streamline Integration through Discontinous Fields.
IEEE-TVCG , 14:680–692, 2008.[28] D. Walfisch, J.K. Ryan, R.M. Kirby, and R. Haimes. One-sided Smoothness-IncreasingAccuracy-Conserving filtering for enhanced streamline integration through discontinuous fields.
Journal of Scientific Computing , 38:164–184, 2009.[29] Eleuterio F. Toro. Riemann Solvers and Numerical Methods for Fluid Dynamics.
Springer-Verlag, Berlin , 2009.[30] A. Harten. High resolution schemes for hyperbolic conservation laws.
Comput. Phys. , 49:357–393, 1983.[31] A. Harten. Adaptive multiresolution schemes for shock computations.
Comput. Phys. , 115:319–338, 1994.[32] B. Costa and W.S. Don. Multi-domain hybrid spectral-WENO methods for hyperbolic conser-vation laws.