For a wide range of partial differential equations, their solutions can be discontinuous. For example, for singularly perturbed problems, when the transient layer is very sharp, the solution can be viewed as discontinuous, see RST:08 ; FR:11 . For linear and nonlinear hyperbolic equations, the discontinuous solutions can be caused by a discontinuous initial data or shock forming, see HW:08 ; Hesthaven:18 . Sometime, the location of discontinuity is known, for example, the discontinuity in the initial or boundary data. For many other cases, the exact location of discontinuity is unknown. When different numerical methods are used to solve such problems, the numerical solutions are often oscillatory near a discontinuity, i.e., it overshoots (and undershoots). In spectral methods and Fourier analysis, it is also called the Gibbs phenomena, see HH:79 ; GS:97 . Such spurious oscillations are unacceptable for many situations, see for example, page 209 of Hesthaven Hesthaven:18 . Such oscillations are essentially caused by a simple fact: continuous functions are used to approximate a discontinuous function. In numerical methods, some tricks are used to eliminate or reduce the Gibbs phenomena: artificial viscosity, limiters, filters, and ENO/WENO schemes, see HW:08 ; Hesthaven:18 ; GS:97 .
On the other hand, it is well known that adaptive finite element method with mesh refinements based on a posteriori error estimation is a very useful procedure to detect the regions and locations with bad approximations, seeVer:13 . A common adaptive finite element algorithm contains the loops of the following steps:
The SOLVE step solves the numerical problems. The ESTIMATE step computes the a posteriori error estimates and indicators. The MARK step marks these elements with big error indicators or small error indicators. The REFINE step refines those elements marked with big error indicators. For some advanced algorithms, an extra COARSEN step is applied to combine these elements with small error indicators. If the solution is discontinuous and often badly approximated, an accurate a posteriori error estimator can often mark these badly approximated elements..
But can adaptive finite element methods get accurate approximations of discontinuous solutions of PDEs without overshooting? For this question, there are some confusing and conflicting folklores: 1. an adaptive mesh cannot reduce the overshootings. This is the common opinion in the numerical hyperbolic equations community with finite difference/finite volume methods, see BO:84 ; Shu:98
. 2. When the mesh is fine enough, the discontinuity is resolved, the overshooting can be reduced. This opinion is probably more common in the adaptive finite element method community wherethe adaptivity can almost cure everything is believed. 3. Since the solution is discontinuous, if methods based on discontinuous piecewise polynomials, for example, discontinuous Galerkin (DG) methods, are used, the overshootings can be eliminated or reduced. For example, in a numerical computation for a singularly perturbed problem in a SIAM review paper FR:11 (p. 165), the authors claimed the DG methods will have no overshootings for almost discontinuous solutions.
In this paper, we want to clarify this question by a simple model problem. Since most of the finite element or DG methods are based on projections or pseudo-projections, the best result a numerical method can achieve of a certain approximation often cannot be better than its -projection, i.e., if is the discrete approximation space, the best we can hope for the numerical error is often in the form (many methods cannot have such a good result, for example, the DG method for the linear hyperbolic equation):
where is the projection operator onto the discrete approximation space . For a step function , we try to study the properties of under a refined mesh to answer the above question.
Based on the computations, we find the following results. For the simple discontinuity-aligned mesh case, piecewise discontinuous approximations are always good. For the general non-matched case, we explain that the piecewise discontinuous constant approximation combined with adaptive mesh refinements is the best choice to achieve accuracy without overshooting. For discontinuous piecewise linear approximations, non-trivial overshootings will be observed unless the mesh is matched with discontinuity. For continuous piecewise linear approximations, non-trivial overshootings will always be observed under regular meshes. We calculate the explicit overshooting values for several typical cases.
We also find that a special high ratio mesh obtained by adaptive mesh refinements combined with coarsening can also reduce the overshooting. But such case is not always possible in two and three dimensions or in one dimension with other error sources.
The paper is organized as follows. Section 2 describes a model problem of approximating a step function by projections onto piecewise constant and linear function spaces. Discontinuous piecewise constant and linear approximations are explicitly computed for the model problem in Section 3. In Section 4, we compute continuous piecewise linear approximations for the model problem under a ”far away assumption”. In Section 5, linear conforming, linear DG, and lowest order mixed finite element methods are tested on a one dimensional singularly perturbed problem. Seveal test problems of 2D linear transport problems with P0- and P1-DGFEMs are done in Sections 6. In Section 7, we make some concluding remarks.
2 A Model Problem
Consider a step function :
Here, the discontinuity gap is . For more general cases with a gap , , it is easy to see the corresponding results are also proportional to .
Let be a one-dimensional mesh on with elements (intervals in the 1D setting) denoted by . We suppose and are big enough positive numbers with respect to the mesh size. Let be the space with polynomials defined on whose degree is less or equal to integer . We consider three approximation spaces defined on :
Let be , , or , the projection is defined by:
We use notations , , and to denote the -projections of on , , and , respectively.
3 Discontinuous Approximations
In this section, we consider two discontinuous approximations: and approximations. Note that, for a discontinuous projection, the approximation is exact on those elements which are away from the discontinuity. Thus, we only need to discuss the special interval that contains the discontinuity.
Let and , we will consider the interval , see (a) of Fig. 1.
3.1 Discontinuous piecewise constant approximation
A simple computation shows that the constant projection of on the interval is:
It is obvious that for , thus there is no overshooting. If or , which corresponds to the cases that the discontinuity is matched with the mesh, the numerical approximation is exact.
Compute the error in the interval , we have,
The worst case is , when . This also matches the a priori analysis that for for an arbitrary small .
With a good a posteriori error estimator that can identify the bad approximated elements, the discontinuity-crossing elements will be found and divided. The and other integration based norms of error will become smaller and smaller. Thus, for this case, the adaptive finite element methods do can achieve accuracy without overshooting.
3.2 Discontinuous piecewise linear approximation
In this subsection, we consider the projection of onto a linear function on . To this end, let and be the linear Lagrange basis functions define on with , , where or , and and .
Let the projection . The projection problem is then
Note that all terms of the matrix problem can be computed exactly. Solve the projection problem, we get
It is easy to show that
The overshooting value can be defined as the following function:
We plot the value of for on the right of Fig. 1. It is easy to see that only when the mesh is aligned with the discontinuity ( or ), is zero. For away from or , the overshooting phenomena is severe. The maximum overshooting value appears at or .
We also notice that if the bisection is used, the relative position of the discontinuity normally will not converges to the left or the right point of the internal. For example, consider the the interval , i.e., the discontinuity is at the position of the internal, if we bisect the interval, the new interval contains the discontinuity is , the discontinuity is at the position. Keeping doing the bisection, we will find that the discontinuity will jump between the and positions. In this case, will always be or , the overshooting value will not decrease. For other initial positions other than the matched case, the overshooting values will oscillate between , for some .
Compute the error in the interval , we have,
When or , the error is zero. The largest error appears at .
Note that for , thus we have
Thus, for discontinuous approximations, even though that the -norm of error of the approximation is better than the approximation, they are of the same order if the mesh is not aligned, and the overshooting cannot be avoided in most cases.
4 Continuous Piecewise Linear Approximation
For the continuous piecewise linear approximation, the exact solution is not in the approximation space even for the discontinuity matched case, and the projection is not local to one element anymore. Luckily, based on the numerical experiments, we can safely assume the following far away assumption for the continuous piecewise linear approximation.
Approximation on nodes far away from the discontinuity (far away assumption): For a node of the mesh , if there is at least one node between the discontinuity position and , the value of at is very close to the exact value of .
4.1 Discontinuity matched mesh
We first consider a possible ideal case: the mesh is aligned with the discontinuity and the mesh is symmetric with respect to the discontinuity at .
Let be a symmetric mesh on with nodes: , where and .
We first test on a uniform mesh with and . We observe a consistent pattern is observed fore different .
From the numerical tests, for example, and in Fig. 2, we note that the values of at the center nodes are identical for different choice for . For big enough, is almost identical to . The values for are , , , , , , , , and . The severe overshootings only appear on the two nodes around the discontinuity. The far away assumption is clearly true in this case.
With the far away assumption, a simple calculation can be done to show why the values of at are close to . We simplify the calculation by assuming that at all , . We assume the values of , where is a small number compared to . We have the rest three values at and to be computed by the projection. By the symmetry of nodes, it is easy to see that at is and . Thus we only need to determine the value of . Let and be the linear Lagrange basis functions on the mesh, with , . Then on can be written as
The value of can be obtained by a simple projection:
A simple calculation shows
Thus, the overshooting value is always close to of the discontinuity jump (2 in our example). The exact value is due to the fact is about in the example.
The errors in the two adjacent elements of is also easy to compute,
4.2 Discontinuity non-matched mesh
In this subsection, we consider the case that that the mesh is not matched with the discontinuity. From the above discussion, we notice that the size of each element is not essential (their ratio is more important, see the example below).
4.2.1 A local uniform mesh
Suppose and , we consider the following mesh with 5 elements with size : , , , , , and , see the left of Fig. 3 for a test mesh with and . The discontinuity cuts through the central interval at . Let be the linear Lagrange basis function for . As before, we assume that is exact at (we omit the small perturbation for simplicity). Then on the interval ,
Solve the following projection problem,
The overshooting value is defined as the following function:
We plot the overshooting value with respect to on the right of Fig.3. The overshooting values lies in the interval
Also, we find that when , the maximum overshooting happen at and when , the maximum overshooting happens at , that is, the maximum overshooting appears at the one of the endpoints of which is farther from the discontinuity. In any choice of , the overshooting is non-trivial.
We have the following error:
For , we get is between and . Note this error matches with the result of (4.3), which corresponds to or .
4.2.2 A local adaptive mesh
We also test the following mesh: , , , , , and , with sizes , , , and , respectively. This kind mesh can be viewed as an example that the mesh is obtained by adaptive refinements.
With this mesh setting and assuming the values of at are exact, by a similar procedure, we get
We show the result on the right of Fig.4, the overshooting value lies in the interval
Also, as before, we find that the maximum overshooting happens at when , and it happens at when . For any choice of , the overshooting is non-trivial. We can vary the element size of elements adjacent from , , , to , and get similar results of non-trivial overshootings. This shows that adaptive mesh refinements by bisection will not make the overshooting phenomena disappear for the continuous linear approximation.
We have the following error:
For , we get is between and .
4.3 Adaptive mesh with coarsening
For the one dimensional problem, there is one special case that the overshooting phenomena can be eliminated.
Consider the following mesh , with , , and , with and . Similar as the uniform grid case, the projection on can be written as
Computed by projection, we get
This means when the ratio is big enough, the overshooting value can be reduced to a very small value. For a one dimensional mesh, this is possible by combining the refining and coarsening, i.e., when an a posteriori error indicator on an element is big, then element is refined, while for those elements with small indicators, they are combined with the close elements to form bigger elements.
We test the case with the following example: the domain is , the initial mesh is . We use the exact error as the error estimator, i.e., , and set the refinement/coarsening criteria as ”refine those elements whose and coarsen those elements whose ”. For this very simple problem, the final mesh will only contains three elements, a center very small element around , and two big elements on the left and right, receptively. After 14 iterations, the final mesh is , and the numerical solution at four nodes are . The overshooting phenomena is almost invisible.
Note that such case will only happen in one dimensional setting. In two and three dimensions, if we have to keep the mesh conforming and regular, the radio of mesh sizes between two adjacent elements sharing a same edge/face cannot be big, thus there will always be non-trivial overshootings. Even for one dimensional problems, such high ratio mesh is probably not possible due to other approximation errors to keep the mesh from coarsening, see a numerical test in Section 5.
5 Numerical Test 1: a 1D singularly perturbed problem
Consider the following problem:
We choose . For this small, the solution is almost identical to , and has an extremely sharp layer at . For the case that mesh size is not small enough to match with the layer, can be viewed as a discontinuity location.
We test the problem with three numerical methods: the linear -conforming finite element method, linear discontinuous Galerkin finite element method, and the lowest-order mixed method.
Let be a one-dimensional mesh of , with and .
Linear conforming finite element method (P1 conforming): Find , such that
where . Note that for a very small , this is the -projection to we discussed.
Linear discontinuous Galerkin finite element problem (P1-DG): Find , such that
where is a big enough number (we can safely choose ). This is the -projection to for a very small .
To introduce the mixed method, we let , then . The mixed variational formulation is: find , s.t.,
We choose the approximation spaces to be and for and , respectively.
Mixed finite element problem (P0 mixed): Find , s.t.,
This is the -projection to for a very small .
In the computation, we use the robust error estimator for the conforming finite element method in Ver:98 ; Kun:02 to drive the adaptive mesh refinement and compute the conforming, DG, and mixed finite element solutions on the same adaptive mesh.
Discontinuity matched mesh case. An initial mesh is chosen. We refine it for 20 times with the maximum refinement strategy . The center nodes are . Numerical approximations with linear conforming, P0-mixed, and P1-DG approximations are shown in Fig. 6. The overshooting value for the linear conforming finite element approximation is , which matches discontinuity matched case. The P0-mixed solution has no overshooting due to that it is a piecewise constant approximation, and the P1-DGFEM has no overshooting since the discontinuity is matched with the mesh.
Discontinuity non-matched mesh case. Next, we choose an initial mesh . With bisections, the discontinuity location will never be exactly matched. The center nodes are . The values of linear conforming solution at two nodes adjacent to are and . This also matches the discussion in (4.4). A value of is found of P1-DGFEM approximation at the node . The value is a little big larger than the discussion in Section 3.2, where piecewise linear function is used to approximate a discontinuous step function, while here it is a discontinuous piecewise linear function. The discontinuous mixed solution has no overshooting for this case.
For a similar problem, in the bottom page 165 of FR:11 , the authors claimed that the DG method for this problem will have no oscillations. This is incorrect, the DG method will have no oscillations unless the discontinuity is matched with the mesh. The numerical test in FR:11 uses a uniform mesh with , which is exactly the discontinuity matched case. Thus, the non-overshooting in FR:11 is not only due to the use of DGFEM, but because of using DGFEM on a discontinuity-matched mesh.
Adaptive mesh with both refinement and coarsening. Next we test the problem with both refine and coarsening. The initial mesh is chosen to be a uniform mesh with , and the nodes on the final mesh are . The conforming finite element solution at these nodes are with almost no overshooting, see Fig. 7.
We should point out that it is not always possible to get a high ratio mesh such as that in of Fig. 7 with a combination of an adaptive refinement and coarsening. We test the problem with the righthand side:
Note that for the linear finite element approximation, there is some approximation error on all elements. Thus these elements are not coarsened, we eventually have a standard adaptively refined mesh and overshooting phenomena appears, see Fig. 8.
6 Numerical Test 2: DG methods for a linear transport equation in 2D
Consider the following linear transport (advection) equation:
where the advective velocity field
is a vector-valued function defined onand . The inflow part of , , is defined in the usual fashion, where denotes the unit outward normal vector to at .
The DGFEM for the transport equation is BMS:04 : find , or such that
where the bilinear form and linear form are defined by
where the term is the usual upwind flux, meaning that takes the value on the inflow side of the edge/face . The sets , , and are the sets of inflow, interior, and outflow edges, respectively. We call the methods P0-DGFEM and P1-DGFEM for and , respectively.
We use a residual type of error estimator similar to that developed in Burman:09 to drive the adaptive mesh refinement.
6.1 A discontinuous solution on a non-matching mesh
Consider the problem: with . The inflow boundary is . Let and . Choose the inflow boundary condition such that the exact solution is
The initial mesh is shown on the left of Fig. 9. The point and are the bottom central node and the top central node, respectively. The mesh and its subsequent bisection mesh are not aligned with the discontinuity although the inflow boundary condition is matched.
On the center of Fig. 9, a final adaptive mesh is shown with many refinements along the discontinuity since the mesh is not matched. On the right of Fig. 9, we show the overshooting values of the P1-DGFEM, with the overshooting value defined as:
From the figure, we see that the overshooting values are oscillate between 0.15 and 0.35 and never goes to zero, which is matched with our analysis. For the P0-DGFEM, the overshooting value is in the order of machine accuracy, , so there is no overshooting at all.
Since the exact solution of the problem is not changing with respect to the -coordinate, we plot a projected solution by plotting the numerical solutions the element center ()/ nodes () suppressing the -coordinate. On Fig. 10, we show the P0-DGFEM projected solution computed on the final adaptive mesh on the left and P1-DGFEM one on the right. We clear see the overshooting of P1-DGFEM.
6.2 Curved transport problem 1
Consider the problem on the half disk . The inflow boundary is . Let , with is the polar angle. Let , , and the inflow condition and the exact solution be
The initial mesh is shown on the left of Fig. 11. On the center of Fig. 11, a final adaptive mesh is shown. On the right of Fig. 9, we show the overshooting values of the P1-DGFEM, with the overshooting value defined as:
From the figure, we see that the overshooting values are oscillate between 0.3 and 0.6 and never goes to zero with the mesh refinements.
We plot projected solutions by plotting the numerical solutions the element center ()/ nodes () with respect to the radius on Fig. 12, we show the P0-DGFEM projected solution computed on the final adaptive mesh on the left and the P1-DGFEM one on the right. We clearly see the overshooting of P1-DG approximation but non-overshooting of the P0-DG approximation.
6.3 Curved transport problem 2
Consider the following problem: with , , and . The inflow boundary is , i.e., the west and north boundaries of the domain. Choose such that the exact solution is
When , the layer is never fully resolved in our experiments and can be viewed as discontinuous. In Fig. 13, we show numerical solutions by P0- and P1-DGFEMs. It is easy to see that there is no overshooting for the P0-DGFEM solution, and non-trivial overshooting can be found for the P1-DGFEM solution.
Similar results are available for these three numerical tests with flux based least-squares finite element formulations can be found in LZ:18 .
7 Concluding Remarks
In this paper, we studies the behavior of using adaptive continuous and discontinuous finite elements to approximate discontinuous and nearly discontinuous PDE solutions.
For a singularly perturbed problem with a small diffusion or other equations with a transient layer, when the mesh is fine enough around the transient layer, the numerical solution can accurately approximate the PDE solution. This is the case that the belief that the adaptive finite element works come from.
For a pre-asymptotic mesh of such layer problems, the solution can be viewed as discontinuous and we can treat it identically as the genuinely discontinuous case.
To approximate a discontinuous solution, the global continuous function is always a bad choice since it will always have non-trivial oscillations. Other than that, it will hard to impose a discontinuous boundary boundary condition with a continuous approximation space. If the exact location of the discontinuity is known, then high order discontinuous approximations on a discontinuity-matched mesh is a good choice. For singularly perturbed problems with a small diffusion, this mesh-matched case is reduced to that the mesh interface lies within the small transient layer. If the location of the discontinuity is unknown, the only way to avoid the overshooting is to use piecewise constant approximations for those elements where a discontinuity cuts through. Discontinuous piecewise linear function approximations will have non-trivial overshooting for most cases. The adaptive finite element method with a good a posteriori error estimator can reduce the -norms (or and other integration based norms) of the error, but it cannot reduce the non-trivial oscillations unless a piecewise constant approximation is used in the discontinuity-crossing elements.
The non-trivial overshooting we find in the paper is of the same magnitude as the Gibbs phenomena in the Fourier analysis, where a number around of the jump gap is found, see HH:79 .
With the above results in mind, when we design finite element methods for problems with unknown location discontinuous solutions, an hp adaptivity approach is preferred. We use high order finite elements (continuous or discontinuous) in the smooth region and use piecewise constant approximations in the discontinuity-crossing elements. This is one of our ongoing work. It is also worth to mention that the method proposed by Lax Lax:06 is exactly such a method in the context of the finite difference method.
- (1) M. J. Berger and J. Oliger, Adaptive mesh refinement for hyperbolic partial differential equations, Journal of Computational Physics, Volume 53, Issue 3, March 1984, Pages 484–512.
- (2) F. Brezzi, L. D. Marini, and E. Süli, Discontinuous galerkin methods for first-order hyperbolic problems, Mathematical Models and Methods in Applied Sciences, 14 (2004), pp. 1893–1903.
- (3) E. Burman, A posteriori error estimation for interior penalty finite element approximations of the advection-reaction equation, SIAM J. Numer. Anal, Vol. 47, No. 5, pp. 3584–3607.
- (4) S. Franz and H-G. Roos, The capriciousness of numerical methods for singular perturbations, SIAM Review, Vol. 53, No. 1 (2011), pp. 157–173.
- (5) D. Gottlieb and C.-W. Shu, The Gibbs phenomenon and its resolution, SIAM Review, 39 644-668, 1997.
- (6) J. S. Hesthaven, Numerical Methods for Conservation Laws: From Analysis to Algorithms, SIAM, 2018.
- (7) J. S. Hesthaven and T. Warburton, Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications, Texts in Applied Mathematics (TAM), vol. 54, Springer, 2008.
- (8) E. Hewitt and R. Hewitt, The Gibbs-Wilbraham phenomenon: An episode in Fourier analysis, Archive for History of Exact Sciences. 21 (2): 129–160, 1979.
- (9) G. Kunert, A note on the energy norm for a singularly perturbed model problem, Computing, vol 69, 3, 2002, 265–272.
- (10) P. D. Lax, Gibbs phenomena, J. Sci. Computing, 28(2-3):445–449, 2006.
- (11) Q. Liu and S. Zhang, Adaptive least-squares finite element methods for linear transport equations based on an H(div) flux reformulation, arxiv 1807.01524, 2018.
- (12) H.-G. Roos, M. Stynes, and L. Tobiska, Robust Numerical Methods for Singularly Perturbed Differential Equations, Springer Ser. Comput. Math., Springer, Berlin, Heidelberg, 2008.
- (13) C.-W. Shu, Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws, in Advanced Numerical Approximation of Nonlinear Hyperbolic Equations, B. Cockburn, C. Johnson, C.-W. Shu, E. Tadmor, and A. Quarteroni, eds., Lecture Notes in Math. 1697, Springer-Verlag, Berlin, 1998, pp. 325–432.
- (14) R. Verfürth, Robust a posteriori error estimators for a singularly perturbed reaction-diffusion equation, Numer. Math. (1998) 78, 479–493.
- (15) R. Verfürth, A Posteriori Error Estimation Techniques for Finite Element Methods, Oxford University Press, Oxford, UK, 2013