Stable and Accurate Filtering Procedures

High frequency errors are always present in numerical simulations since no difference stencil is accurate in the vicinity of the π\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\pi $$\end{document}-mode. To remove the defective high wave number information from the solution, artificial dissipation operators or filter operators may be applied. Since stability is our main concern, we are interested in schemes on summation-by-parts (SBP) form with weak imposition of boundary conditions. Artificial dissipation operators preserving the accuracy and energy stability of SBP schemes are available. However, for filtering procedures it was recently shown that stability problems may occur, even for originally energy stable (in the absence of filtering) SBP based schemes. More precisely, it was shown that even the sharpest possible energy bound becomes very weak as the number of filtrations grow. This suggest that successful filtering include a delicate balance between the need to remove high frequency oscillations (filter often) and the need to avoid possible growth (filter seldom). We will discuss this problem and propose a remedy.


Introduction
For reliable solutions to initial boundary value problems (IBVP), stability is required. This can be achieved by using schemes on SBP form together with weak imposition of boundary conditions. A central feature with these schemes is that the discrete spatial operator is associated with a corresponding integration procedure (inner product). Stable and high order accurate discretizations of SBP type have been known for a long time [10].
High frequency errors are always present in numerical simulations since no difference stencil is accurate in the vicinity of the π-mode. To remove the defective high wave number information from the solution, artificial dissipation operators or filter operators may be B Jan Nordström jan.nordstrom@liu.se Tomas Lundquist tomlu35@gmail.com applied. Artificial dissipation operators preserving the accuracy and energy stability of SBP schemes are available [6,9]. However, for filtering procedures it was recently shown that stability problems may occur, even for originally energy stable (in the absence of filtering) SBP based schemes [7].
The application of an oscillation-reducing filter during a numerical simulation can be seen as one particular type of so-called transmission problems recently studied in [7]. It was shown that in order to preserve stability of the underlying numerical scheme, the transmission problem must be contractive with respect to the same inner product for which the discrete spatial operator is semi-bounded. Unfortunately, the boundary closures for explicit filters available in the previous literature do not satisfy this contractive property even in the most simple low order case. On the other hand, for transmission problems involving adaptive mesh refinement, so-called SBP preserving interpolation operators sometimes satisfy the required contractive property [5]. In this paper we will demonstrate that the SBP preserving concept (which we will denote inner product preserving (IPP) in the rest of the paper) is central also for the problem of constructing stable filters. In particular we will show that an IPP property is necessary but not sufficient for an explicit filter to be contractive. Sufficiency will be indicated by a set of separate criteria for filters of finite difference type. If the same filter operators are implemented in an implicit way, we prove that the IPP property is in itself also sufficient for stability. This paper is organized as follows. The filter transmission problem is introduced in Sect. 2, together with the stability condition for explicit filters derived in [7]. Explicit finite difference stencils and a previous method for constructing boundary closures are reviewed in Sect. 3. A necessary IPP condition for stability is then formulated in Sect. 4, which explains the boundary instabilities previously observed. A solution to the stability problem is proposed in the form of a new filter construction formula with the IPP property. The new stable filters are extended to curved multi-dimensional geometries in Sect. 5. In Sect. 6 we investigate the numerical properties of the new explicit and implicit filter implementations, demonstrating the stability and excellent oscillation-reducing properties for a boundary layer problem. Finally in Sect. 7, we draw conclusions.

The Transmission Problem
In this section we briefly revisit the setting in [7] and introduce the necessary notation.

Continuous Model Problem
Consider a possibly nonlinear partial differential equation governing the evolution of a solution field u(t, x) on the domain Ω ∈ R d . We write, where D is a differential operator, and L is a boundary operator. We assume that the problem is semi-bounded, i.e. Ω φ D(φ) ≥ 0, for smooth functions φ satisfying the homogenous version (g = 0) of the boundary conditions in (1). Next, we semi-discretize the problem with homogeneous boundary conditions, and obtain, where x is a vector containing the grid points. The operator D approximates D and includes the weakly imposed boundary conditions. Further, we assume that the discrete operator is also semi-bounded, i.e. that where P is a diagonal operator (inner product and norm) containing the quadrature weights associated with D. Stability then follows, since, holds.
Next, we consider the problem of applying an oscillation-reducing filter to the solution at some arbitrary point t 1 in time. This operation can be seen as a so-called transmission problem [7]. Building directly on the semi-discretized form of the equations, we consider the problem, where Ψ is a general filter function which can be either implicit or explicit (Ψ = Ψ [U (t 1 )]).
For stability we require the contractive property, which by (3) yields the desired result, The estimate (6) shows that (4) is estimated in the initial data of the problem (1). We say that schemes that satisfy (6) are time-stable. It was shown in [7] that schemes that are not time-stable can lead to erroneous energy growth. We conclude this section by considering the special case of a linear explicit filter, where F is a matrix (see Fig. 1). For such filters, (5) implies that the contractive property of F first identified in [7], must hold. Note that in (7) there is no distinction made between explicit filters on constant coefficient or semi-linear (i.e. F = F (U )) form. In both cases, (7) leads directly to time- Fig. 1 Schematic of semi-discrete filter problem

Explicit Finite Difference Filters
We specify x = (0, Δx, 2Δx, . . . , 1) to be a uniformly spaced grid vector in 1D with N + 1 grid points, i.e. Δx = 1/N , and let D be a high order SBP finite difference operator with the inner product P. The highest solution frequency which can be resolved on the grid is the so-called π-mode, The central feature of most oscillation-reducing explicit filters is to remove the π-mode from a numerical solution through the application of a dissipative high order derivative approximation.
It is well known that, down to a scaling factor, compact stencil approximations of even order derivatives leave the π-mode intact. For example, a symmetric second order accurate stencil approximating the 2n:th order derivative satisfies, see e.g. [2]. Since D 2n is itself scaled with a factor of 1/Δx (2n) , second order accuracy means that polynomials up to order 2n + 1 are evaluated exactly. In particular, polynomials up to order 2n − 1 are then differentiated exactly to zero, where the power j should be interpreted as an element-wise operation on x. It follows that a filter defined by F = I − α(−1) n D 2n , where α = Δx (2n) /2 (2n) , simultaneously cancels the π-mode while leaving low order frequency modes unaltered to within 2n:th order of accuracy, These filter properties only apply in the interior of the domain. In [2], boundary closures were constructed for stencils up to sixteenth order (n = 8), leading to operators D 2n for a bounded domain satisfying the symmetry and accuracy conditions, Based on these operators, an explicit filter defined by F = I − α(−1) n D 2n then satisfies, where in addition (10) and (11) applies to the interior stencil of F . The fact that the spectral radius of F is smaller than unity implies that the contractive property (7) holds if P is proportional to the unit matrix, i.e. if P = Δx I . Unfortunately, this is typically not the inner product for which semi-boundedness (2) can be proven. For other P:s, (7) can not be guaranteed based on the design requirements in (12) and (13). As observed in [7], even in the most simple and lowest order case, the operators in [2] fail to satisfy (7). Consider for example a second order accurate SBP norm with N = 3 and Δx = 1/3, A symmetric filter of the form (14) and (15) for n = 1 is given by, The eigenvalues of F T PF − P become [− 0.9375, − 0.5890, − 0.1250, 0.0265]. The positive eigenvalue indicates that the numerical scheme is not time-stable, which can lead to unwanted energy growth.

Stable Filters: Theory and a New Operator Definition
Even though in this paper we focus on filters to be used in conjunction with finite difference schemes, we stress that the formulation of the transmission problem in Sect. 2 is general and applies to all types of semi-bounded schemes. In this section we will derive stability results for general filters on both explicit and implicit form. We will then apply this theory to propose a new and improved finite difference filter definition.

A Necessary Stability Condition for Explicit Filters
The contractivity condition (7) implies that the construction of filter boundary closures must relate to the inner product P. In Proposition 1 below we will show that the only way to achieve (7) is to let F be associated with another operatorF of the same accuracy, such that they form a so-called IPP pair, defined by, This might seem like a counter-intuitive statement. The operatorF does not appear in the scheme (4), but somewhat surprisingly, its existence is necessary for the transmission problem to be stable. To prove this, we need the following lemma. Proof See Lemma 14 in [8].
We can now prove Proposition 1 Let the operator F in (4) be n:th order accurate, i.e. it satisfies (15). In order for the contractivity condition (7) to hold with respect to the inner product P, it is necessary that also the operatorF in (17) is n:th order accurate.
Proof Assume that (15) holds, and consider the quadratic form, for j = 0, 1, . . . , n − 1. It now follows from Lemma 1 that F T PF − P must be an indefinite matrix unless we also have, Inserting (15) again and then multiplying with P −1 , this yields, i.e.F is also n:th order accurate.
For a filter defined as before by F = I − α(−1) n D 2n , (19) does by no means follow from the properties (12) and (13). In fact one can easily verify numerically that (19) in general does not hold, thus showing that the filter definition F = I − α(−1) n D 2n together with (12) and (13) is inadequate. With Proposition 1 in place, we are now equipped to reconsider the problem of constructing stable filter boundary closures. Indeed, we can rule out all cases whereF as defined in (17) does not satisfy (19).

A Modified Finite Difference Filter Definition
By taking advantage of the fact that high order derivative operators D 2n satisfying (12) and (13) are already available, we consider the new filter definition, where P = ΔxH.

Remark 1
This modification is similar to the one proposed in [6] for artificial dissipation.
Note that the new definition (20) satisfies the same accuracy conditions (15) as before. Indeed, (13) implies that F reduces to an identity operator when applied to polynomials of order n−1.
Moreover, from the symmetry of D 2n we have, i.e. the matrix PF is also symmetric. In the IPP definition (17), this leads tõ i.e. F forms an IPP pair together with itself. The condition (19) necessary for stability hence follows directly by construction from the accuracy of F itself in (15). Note that while (19) is not sufficient in order to prove that contractivity (7) is satisfied, it is straightforward to verify numerically if this is the case or not, for a given combination of F and P.
For example, n = 1 together with the second order accurate norm (16) yields,  (7) is not, for example by adjusting the quadrature weights in (16). To make the theory complete, we have developed a set of sufficiency criteria to test whether a given filter of the type (20) is contractive or not. These criteria are independent of the grid size N , see "Appendix A". Using these criteria we can prove that no eigenvalue of F T PF − P is positive for standard SBP finite difference norms up to eighth order accuracy.
To summarize this section, Proposition 1 was applied to propose a plausible candidate for a new filter definition (20). Furthermore, for each combination of operators P and D 2n , the contractivity condition in (7) can be verified using the criteria for sufficiency given in "Appendix A".

A Provably Stable Implicit Implementation
As we have seen, the IPP property (17) for explicit filters defines a new operatorF which is not present in the original transmission problem (4). However, using an implicit formulation we can use this operator to design a provably contractive scheme. In (4), consider the previous explicit filter with an added implicit correction term, whereF is defined as in (17). Note that if (19) holds, then the accuracy of this condition is of the same order as the operator F itself. As we shall now prove, (17) in combination with (22) is both necessary and sufficient for contractivity. (17), the estimate

Proposition 2 Consider the implicit filter implementation (22). By definingF as in
Proof Multiplying (22) with V (t 1 ) T P from the left yields, , U (t 1 )) P . After adding and subtracting U (t 1 ) 2 P on the right hand side, we finally get, which concludes the proof.
We recall that the contractivity property (5) together with the estimate (3) of U leads to time-stability (6) of the transmission problem.

Extensions
To accommodate for non-Cartesian geometries in multiple dimensions, SBP finite difference discretizations of IBVPs can be extended through a combination of tensor products and curvilinear transformations in such a way that semi-boundedness (2) is preserved. In this section we will discuss how such extensions relate to filters.

Multiple Dimensions
Finite difference formulations in 1D can be extended to multi-dimensional Cartesian grids by the application of tensor products. For example, in the 2D case, let where ⊗ denotes the Kronecker (or tensor) product, and where P x , P y , F x and F y are integration and filter operators in the two coordinate directions, respectively. The IPP property (17) immediately extends from 1D to 2D operators, sincẽ following elementary rules of the Kronecker product. This means that the implicit filter implementation (22) extends directly to multiple dimensions in a stable and accurate way. As for the contractive property (7), consider By adding and subtracting 1 2 Since F T x P x F x + P x and F T y P y F y + P y only have positive eigenvalues, it follows that (7) is satisfied provided that it is satisfied for the 1D operators, i.e.
We have thus shown that stability is preserved for both the explicit and implicit implementations in 2D. By induction, these arguments can be extended to any number of dimensions.

Curvilinear Coordinates
Letx denote a Cartesian grid in any number of space dimensions, and let x = x(x) denote a curvilinear transformation such that we obtain a new grid vector x. Moreover, let J be the diagonal matrix obtained by inserting values from the Jacobian determinant of this transformation at the grid points. From a discrete inner productP defined on the Cartesian grid, we can define a new operator P = JP containing the positive quadrature weights on x.
We define the new filter operator on x as, whereF is constructed with respect to the Cartesian grid exactly as in the previous section. An IPP pair is now obtained by analogously definingF as, showing that the implicit implementation (22) of F in (23) is stable and accurate. Also the contractivity condition (7) extends to the curvilinear case, and we can prove Proof Consider, Since bothP and J are diagonal, they commute, leading to, and the result follows.

Numerical Calculations
We test the numerical properties of the filter definition in (20), henceforth referred to as the new filters, and compare to the previous version with F = I − α(−1) n D 2n , which we refer to as the old filters. In both cases, the high order derivative operators D 2n are taken from [2]. If filtering is carried out repeatedly after each or every few time steps in a simulation, the errors will accumulate over time. To avoid a reduction in convergence rate, the filter accuracy should be at least one order higher than that of the numerical scheme itself. To limit the number of cases considered, we shall for each value of n associate F with a finite difference operator of order 2(n − 1) in the interior, and order n − 1 in the boundary closure. The lower order accuracy in the boundary closure is a theoretical upper limit for the type of operators considered, see [3,4].

Filter Properties
First we investigate the effect of the explicit (4) and implicit (22) filter implementations to functions with various frequencies. Consider the test function, where 0 ≤ ξ ≤ π is the so-called wave number. Note that ξ = π yields the π-mode f π in (8), i.e. the highest frequency that can be resolved on the grid. Consider the uniform grid x in 1D with N + 1 grid points, and define the vector U such that U i = u(x i ). In Fig. 2 we plot U for a few different frequencies before and after filtering. Here we have used N = 32 and n = 2, i.e. a fourth order derivative D 2n together with a norm P associated with an SBP finite difference operator of order 2 in the interior. As expected, filtering results in a cancellation of the π-mode in the interior of the domain, while the damping of low frequencies is small. In order to compare the different methods to one another, it is convenient to first separate the action in the interior of the domain from that of the boundary closures. In Fig. 3 we plot the interior damping as a function of the wave number, using n = 2, 3, 4 and 5. Notice that the new and old explicit filter operators are identical in the interior, which is why no there is distinction between them here. We observe that all filters succeed in canceling the π-mode, and that higher order filters are less dissipative than lower order ones for all wave numbers smaller than π. For the same value of n, the implicit implementation is less dissipative than the explicit one.
Finally, in Fig. 4 we take a closer look at the results from filtering the π-mode close to a domain boundary. The cases considered here do not differ significantly from each other, with the new filters showing slightly better damping for lower values of n, but slightly worse for n = 5.

Boundary Layer Calculation
One possible source for unwanted oscillations in a numerical solution is when a second order derivative is approximated by applying a central first derivative operator twice, leading to a non-compact stencil. See e.g. [1] for an analysis of such oscillatory error modes. Since central difference approximations of odd order derivatives by themselves cancel out the πmode, there is no natural mechanism in such schemes to damp out oscillations once they start appearing.
As a model problem, we consider the advection-diffusion equation with well-posed boundary conditions, An exact steady-state solution with a steep boundary layer at x = 1 is given by, We let = 0.1, and use first derivative finite difference SBP operators to approximate both u x and u x x . More details about this discretization technique can be found in [1]. In Fig. 5 we show examples of error distributions e = U − u e (x), where U is the numerical solution obtained either with or without filtering. We have used finite difference SBP operators of fourth order in the interior (and hence n = 3), and compute the steady-state solution using the classical RK4 method in time. As initial solution we use f(x)=0. The time step size is chosen to be Δt = Δx 2 /4 , and filters are applied after each time step. In Fig. 6 we plot the error in maximum norm (with respect to space) as a function of time, showing that steady-state is reached at around t = 5 in all cases.
The positive effect of filtering is apparent. Dominated by oscillations, the errors are almost eliminated outside the boundary layer. For the new time-stable techniques, the oscillations close to x = 1 are also reduced significantly. The old method on the other hand is showing Convergence of maximum error with and without filtering. The dotted line indicates third order convergence clear problems, and even fails to converge with the same rate as the unfiltered scheme, see Fig. 7 for convergence plots of the error in maximum norm. Since the formal order of accuracy is the same in all cases, we conclude that this convergence issue is due to an instability caused by the old filter.
Since the steepest gradients are found within the boundary layer close to x = 1, it makes sense to reduce the grid spacing in this region. In the final test, we apply the smooth grid transformation, where d = 3/2, andx is a uniform grid with N + 1 grid points. The filter operators are modified according to Sect. 5.2, and results are shown in Figs. 8 and 9. The error levels are significantly reduced compared to the Cartesian grid calculations, and the new filters are again superior to the old ones. Close to x = 1, the errors from the explicit and implicit implementation are comparable. In the interior however, the implicit method is more accurate.

Conclusions
A new inner product preserving (IPP) filter definition for finite difference schemes was proposed. The IPP property was shown to be necessary for the stability of explicit filters in general, while sufficiency was verified for each operator individually. An implicit implementation of the same operators was also proposed, for which the IPP property was shown to be sufficient for stability. These theoretical results extend to filters used in conjunction with other types of schemes as well, and even to nonlinear filters.
The new filters can easily be extended to multiple dimensions and complex geometries in such a way that stability and accuracy is preserved. The advantage compared to a previous unstable filter definition was demonstrated by numerical calculations for an advection-diffusion boundary layer problem, using non-compact stencil second order derivative approximations.
Acknowledgements Open access funding provided by Linköping University.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

A Sufficient Conditions for Contractivity
In this section we derive sufficient conditions for explicit finite difference filters of the form (20) to be contractive (7). Due to symmetry, it is sufficient to consider the problem on a half-line. Thus we denote the diagonal elements in H with h j > 0, j = 0, 1, . . . , ∞. Recall that the inner product H has been scaled such that h j is independent of Δx.
Let D 1 denote the undivided forward difference operator, The symmetric operators D 2n utilized to construct the filters in [2] can all be written on the following form, Thus, we can write (20) as which leads to For contractivity, it is thus sufficient to consider the following matrix associated with each filter, Now, let e j be a vector with a single non-zero value of 1 in the j : th position, and 0 everywhere else, We can now write, The second term on the right hand side of (27) can thus be expanded into the sum of outer products, Next, we will expand the unit matrix in (27) in a similar way, so that each term matches the non-zero structure of the corresponding term (D n 1 e j )(D n 1 e j ) T above. From the definition of e j , it is clear that the vector D n 1 e j is given by the j : th column of the matrix D n 1 , and the structure of D n 1 is illustrated in "Appendix A.1" through A.5. To accomplish the stated goal, we expand I as, where, and, We can now write (27) as, where both terms inside the summation have a non-zero part with the same position and dimension. Moreover, the only eigenvector with a non-zero eigenvalue to the outer product (D n 1 e j )(D n 1 e j ) T is D n 1 e j itself, since all other eigenvectors are given by the orthogonal set to D n 1 e j , The single non-zero eigenvalue λ is thus given by where · denotes the standard euclidean vector norm. Hence, a sufficient condition for M to be negative semi-definite is In terms of the quadrature weights, we have thus derived the following simple test for contractivity.

Proposition 3
The finite difference filter (20) together with (26) is contractive in the sense of (7) if the quadrature weights in H satisfy the following set of inequalities, where D 1 is the undivided forward difference operator in (25), and D n 1 e j is given by the j :th column of D n 1 .
Note that Proposition 3 provides a condition for contractivity which is sufficient but not necessary, i.e. the opposite of the situation in Proposition 1. In other words, a filter might still be contractive even if this condition is not satisfied. For example, in the interior of D n 1 , each column contains the binomial coefficients with respect to n (with positive or negative sign, see the examples in "Appendix A.1" through A.5), and the euclidean norm is given by Vandermonde's identity as, For sufficiently large values of j (i.e. such that both j > n and h j = 1), the condition in Proposition 3 thus reduces to, and by inspection this inequality only holds for values of n between one and ten (i.e. for filters up to order 20). Hence, n = 10 constitutes a hard upper limit above which Proposition 3 is no longer applicable. As we shall see, the condition in Proposition 3 will be sufficient to prove contractivity in the large majority of cases with classical finite difference SBP operators (provided that n ≤ 10). However, there are some cases where it becomes insufficient due to some quadrature weights in the boundary part of H being much smaller than unity. A more powerful test can be obtained by considering the complete boundary part of M as a single matrix, which is given below. Diag(h 0 , h 1 , . . . , h r , 1, 1 . . .) be a given quadrature with r + 1 nonunit weights. Then the associated filter (20) and (26) is contractive in the sense of (7) if n ≤ 10 (for the interior part), and max(r ,n) j=0 − 2 n + 1 I j + 1 h j 2 2n (D n 1 e j )(D n 1 e j ) T ≤ 0,

Proposition 4 Let H =
for the boundary part.
The principal difference between Propositions 3 and 4 is that the latter is less restrictive. The reason for this is that a single eigenvalue problem is considered for the whole boundary contribution in the energy method, whereas in Proposition 3 the contribution is split between each h j individually. The advantage of Proposition 3 is that it is much easier to apply, providing a clearly defined interval for each quadrature weight. In most cases this is enough to verify contractivity. In some cases however, the more powerful test of Proposition 4 is required.

A.3 Sixth Order Case
The sixth order operator (n = 3) uses,

A.4 Eighth Order Case
The eighth order case (n = 4) is given by,

A.5 Tenth Order Case
In the tenth order case (n = 5) we have,

A.6 High Order Quadrature Rules
The non-unit quadrature weights associated with classical finite difference SBP operators of minimal boundary closure dimension are given by, the other hand, Proposition 4 can be applied to prove that contractivity holds even with this quadrature for all filters down to second order (n = 1). For the first three quadratures listed above, we have found that Proposition 3 is sufficient to prove contractivity of all filters of order 20 (i.e. n = 10) or less.