Doubling the convergence rate by pre- and post-processing the finite element approximation for linear wave problems

TL;DR

This paper introduces a novel pre- and post-processing algorithm that doubles the convergence rate of finite element methods for linear wave problems, effective on unstructured meshes and complex domains.

Contribution

A new pre- and post-processing approach that enhances convergence rates without loss on unstructured meshes or near boundaries, requiring minimal additional computational effort.

Findings

Achieves up to 2p order convergence in energy norm.

Effective on unstructured meshes and heterogeneous domains.

Requires solving at most q linear systems for q-step processing.

Abstract

In this paper, a novel pre- and post-processing algorithm is presented that can double the convergence rate of finite element approximations for linear wave problems. In particular, it is shown that a $q$-step pre- and post-processing algorithm can improve the convergence rate of the finite element approximation from order $p+1$ to order $p+1+q$ in the $L^2$-norm and from order $p$ to order $p+q$ in the energy norm, in both cases up to a maximum of order $2p$, with $p$ the polynomial degree of the finite element space. The $q$-step pre- and post-processing algorithms only need to be applied once and require solving at most $q$ linear systems of equations. The biggest advantage of the proposed method compared to other post-processing methods is that it does not suffer from convergence rate loss when using unstructured meshes. Other advantages are that this new pre- and post-processing method is straightforward to implement, incorporates boundary conditions naturally, and does not lose accuracy near boundaries or strong inhomogeneities in the domain. Numerical examples illustrate the improved accuracy and higher convergence rates when using this method. In particular, they confirm that $2p$-order convergence rates in the energy norm are obtained, even when using unstructured meshes or when solving problems involving heterogeneous domains and curved boundaries.

Figures

Tables

Table 1. Table 1 . Number of processing steps α 𝛼 \alpha and β 𝛽 \beta at initial and final time to improve the convergence rate by q 𝑞 q orders.

$q$	1	2	3	4	5	6	7	8	9	10
$\alpha$	1	1	2	2	3	3	4	4	5	5
$\beta$	0	1	1	2	2	3	3	4	4	5

Table 2. Table 2 . Results for the 1D problem showing relative errors in the energy norm of the degree- p 𝑝 p spectral element method with ( q = p 𝑞 𝑝 q=p ) and without ( q = 0 𝑞 0 q=0 ) order- q 𝑞 q improving pre- and post-processing and using N 𝑁 N elements per wavelength.

		$q=0$			$q=p$
	N	$e_E$	ratio	order	$e_E$	ratio	order
	20	1.66e-01			1.54e-01
$p=1$	40	4.90e-02	3.39	1.76	3.86e-02	4.00	2.00
	80	1.71e-02	2.86	1.52	9.65e-03	4.00	2.00
	5	7.23e-02			6.02e-02
$p=2$	10	9.68e-03	7.47	2.90	3.63e-03	16.6	4.05
	20	2.00e-03	4.85	2.28	2.25e-04	16.1	4.01
	5	3.56e-03			4.12e-04
$p=3$	10	4.19e-04	8.50	3.09	6.41e-06	64.2	6.01
	20	5.04e-05	8.31	3.05	1.00e-07	64.0	6.00

Table 3. Table 3 . Results for the 1D problem showing relative errors in the L 2 superscript 𝐿 2 L^{2} - and energy norm of the degree-3 spectral element method with order- q 𝑞 q improving pre- and post-processing and using N 𝑁 N elements per wavelength.

		$p=3$
	N	$e_0$	ratio	order	$e_E$	ratio	order
	5	4.14e-04			6.58e-04
$q=1$	10	6.42e-06	64.4	6.01	3.12e-05	21.1	4.40
	20	1.00e-07	64.0	6.00	1.70e-06	18.4	4.20
	5	4.12e-04			4.25e-04
$q=2$	10	6.43e-06	64.1	6.00	7.35e-06	57.7	5.85
	20	1.00e-07	64.1	6.00	1.41e-07	52.1	5.70

Table 4. Table 4 . Results for the square-domain problem showing relative errors in the energy norm of the degree- p 𝑝 p mass-lumped triangular element method with ( q = p 𝑞 𝑝 q=p ) and without ( q = 0 𝑞 0 q=0 ) order- q 𝑞 q improving pre- and post-processing.

		$q=0$			$q=p$
	mesh	$e_E$	ratio	order	$e_E$	ratio	order
	3	4.02e-02			4.24e-03
$p=1$	4	1.92e-02	2.09	1.07	1.09e-03	3.88	1.96
	5	9.33e-03	2.06	1.04	2.70e-04	4.05	2.02
	1	1.52e-02			4.64e-04
$p=2$	2	3.80e-03	4.01	2.01	2.89e-05	16.1	4.01
	3	9.46e-04	4.01	2.00	1.82e-06	15.9	3.99
	1	1.00e-03			1.16e-06
$p=3$	2	1.26e-04	7.92	2.99	1.75e-08	66.1	6.05
	3	1.58e-05	8.00	3.00	2.97e-10	59.0	5.88

Table 5. Table 5 . Results for the square-domain problem showing relative errors in the L 2 superscript 𝐿 2 L^{2} - and energy norm of the degree-3 mass-lumped triangular element method with order- q 𝑞 q improving pre- and post-processing.

		$p=3$
	mesh	$e_0$	ratio	order	$e_E$	ratio	order
	1	1.68e-06			3.24e-05
$q=1$	2	4.06e-08	41.4	5.37	1.99e-06	16.3	4.03
	3	1.11e-09	36.5	5.19	1.22e-07	16.3	4.02
	1	1.27e-06			5.09e-06
$q=2$	2	1.92e-08	65.8	6.04	1.51e-07	33.7	5.07
	3	3.23e-10	59.6	5.90	4.56e-09	33.2	5.05

Table 6. Table 6 . Results for the circle-domain problem showing relative errors in the energy norm of the degree- p 𝑝 p mass-lumped triangular element method with ( q = p 𝑞 𝑝 q=p ) and without ( q = 0 𝑞 0 q=0 ) order- q 𝑞 q improving pre- and post-processing.

		$q=0$			$q=p$
	mesh	$e_E$	ratio	order	$e_E$	ratio	order
	3	4.21e-01			4.20e-01
$p=1$	4	1.09e-01	3.85	1.95	1.00e-01	4.20	2.07
	5	3.28e-02	3.32	1.73	2.45e-02	4.07	2.03
	1	4.40e-02			2.28e-02
$p=2$	2	7.90e-03	5.57	2.48	1.27e-03	17.9	4.16
	3	1.87e-03	4.23	2.08	7.35e-05	17.3	4.11
	1	2.68e-03			1.75e-04
$p=3$	2	3.35e-04	8.01	3.00	2.22e-06	78.8	6.30
	3	4.20e-05	7.98	3.00	2.96e-08	75.0	6.23

Table 7. Table 7 . Results for the square-domain problem showing relative errors in the energy norm e E subscript 𝑒 𝐸 e_{E} of the degree- p 𝑝 p mass-lumped triangular element method with order- p 𝑝 p improving pre- and post-processing using the conjugate gradient method with N i ​ t subscript 𝑁 𝑖 𝑡 N_{it} iterations. The number of time steps N T subscript 𝑁 𝑇 N_{T} is listed in the third column.

	mesh	$N_T$	$N_{it}=0$	10	100	1000	$\infty$
	3	6190	4.02e-02	5.79e-03	4.22e-03	4.24e-03	4.24e-03
$p=1$	4	12342	1.92e-02	1.67e-03	1.09e-03	1.09e-03	1.09e-03
	5	24727	9.33e-03	4.59e-04	2.85e-04	2.70e-04	2.70e-04
	1	3586	1.52e-02	6.52e-04	4.62e-04	4.64e-04	4.64e-04
$p=2$	2	7222	3.80e-03	9.36e-05	2.91e-05	2.89e-05	2.89e-05
	3	14529	9.46e-04	1.40e-05	1.98e-06	1.82e-06	1.82e-06
	1	7399	1.00e-03	1.19e-04	2.91e-06	1.16e-06	1.16e-06
$p=3$	2	14879	1.25e-04	1.16e-05	4.55e-07	1.77e-08	1.75e-08
	3	29827	1.58e-05	1.38e-06	3.92e-08	7.28e-10	2.97e-10

Equations

$${\partial }_t^{2}u$$

$$u$$

$$u{|}_{t=0}$$

$${\partial }_{t}u{|}_{t=0}$$

$$c_0\le c\le c_1\text{and}{\rho }_0\le \rho \le {\rho }_1$$

$$⟨{\partial }_t(\rho {\partial }_{t}u),w⟩+a(u,w)$$

$$(u,w{)}_{\rho }$$

$$a(u,w)$$

$$U_h=U({\mathcal{T}}_h,\hat{U})$$

$$({\partial }_t^2{u}_h,w{)}_{\rho }+a(u_h,w)$$

$$({\Pi }_{h}u,w{)}_{\rho }=(u,w{)}_{\rho }$$

$$\mathcal{L}:=-{\rho }^{-1}\nabla \cdot c\nabla ,$$

$$a({\mathcal{L}}^{-1}f,w)$$

$$-{\rho }^{-1}\nabla \cdot c\nabla u$$

$$u$$

$$({\mathcal{L}}_h{u}_h,w{)}_{\rho }$$

$$a({\mathcal{L}}_h^{-1}{f}_h,w)$$

$${\partial }_t^{2}u+\mathcal{L}u$$

$${\partial }_t^2{u}_h+{\mathcal{L}}_h{u}_h$$

$${\partial }_t^{k+2}u+\mathcal{L}{\partial }_t^{k}u$$

$${\partial }_t^{k+2}{u}_h+{\mathcal{L}}_h{\partial }_t^k{u}_h$$

$${\partial }_t^{k+2}u$$

$${\partial }_t^{k}u$$

$${\partial }_t^{k+2}{u}_h$$

$${\partial }_t^k{u}_h$$

$$r_k^{(1)}(w)$$

$$r_{k,h}^{(1)}(w)$$

$${\partial }_t^{k+2\alpha }u$$

$${\partial }_t^{k}u$$

$${\partial }_t^{k+2\alpha }{u}_h$$

$${\partial }_t^k{u}_h$$

$$\displaystyle\big{(}r_{k}^{(\alpha)}\circ r_{k}^{(-\alpha)}\big{)}(w)$$

$$\displaystyle\big{(}r_{k,h}^{(\alpha)}\circ r_{k,h}^{(-\alpha)}\big{)}(w)$$

$$u_{0,h}$$

$$v_{0,h}$$

$$u_{h\ast }(T)$$

Full text

The biggest advantage of the proposed method compared to other post-processing methods is that it does not suffer from convergence rate loss when using unstructured meshes. Other advantages are that this new pre- and post-processing method is straightforward to implement, incorporates boundary conditions naturally, and does not lose accuracy near boundaries or strong inhomogeneities in the domain. Numerical examples illustrate the improved accuracy and higher convergence rates when using this method. In particular, they confirm that $2p$-order convergence rates in the energy norm are obtained, even when using unstructured meshes or when solving problems involving heterogeneous materials and curved boundaries.

*The author has been funded by the Austrian Science Fund (FWF) through the project F 65 ‘Taming Complexity in Partial Differential Systems

1. Introduction

When solving wave propagation problems, the finite element method offers a good alternative to the popular finite difference method when the effects of the geometry, e.g. the geometry of objects in scattering problems or the topography in seismic models, need to be accurately modelled. The error of the finite element solution is at most of order $h^{p}$ in the energy norm and order $h^{p+1}$ in the $L^{2}$-norm, where $h$ denotes the mesh resolution and $p$ the polynomial degree of the finite element space. We show that, by carefully discretising the initial values and by post-processing only the final solution, we can improve the convergence rates in both norms up to order $2p$ .

Post-processing methods for the finite element method have been known for several decades and have been applied to numerous problems, including elliptic problems [3, 29], parabolic problems [30], and first-order hyperbolic problems [1, 5]. For an historic overview, see, for example, [5] and the references therein. These methods exploit the fact that the numerical error converges faster in a negative-order Sobolev norm than in the $L^{2}$-norm. Typically, the post-processed solution is obtained by convolving the finite element solution with a suitable kernel. The error of the post-processed solution can then be bounded by the numerical error and difference quotients of the numerical error in a negative-order Sobolev norm.

In this paper, we show that the numerical error for the linear wave equation converges with rate $\min(p+1+m,2p)$ in the $(-m)$th-order Sobolev norm. Difference quotients of the numerical error converge with the same rate in case translation-invariant meshes are used. In case of unstructured meshes, however, the super-convergence rate of the difference quotients, and therefore the post-processed solution, is (partially) lost [9, 21, 17]. We therefore present an alternative post-processing algorithm that fully preserves the super-convergence rate when using unstructured meshes.

The accuracy of this new approach relies on (adapted) negative-order Sobolev norms of time derivatives of the numerical error instead of difference quotients. These time derivatives also play an important role in the error analysis of finite element methods for nonlinear problems [20]. They maintain the full super-convergence rate when using unstructured meshes if the initial values are carefully discretised. To improve the convergence rate by $q$ orders up to a maximum of order $2p$, the proposed algorithm requires solving at most $q$ elliptic problems for discretising the initial value and $q$ elliptic problems for post-processing the final solution. These elliptic problems are solved with a finite element method using the same mesh as for the time-stepping. For post-processing, elements of degree $p+q$ are used instead of degree $p$. The resulting linear systems of equations can be solved with a direct solver or with an iterative solver and for the latter, good initial guesses can be obtained from the unprocessed initial and final values.

The biggest advantage of this new method is that the super-convergence rates are fully maintained when using unstructured meshes. Other advantages are that the method is easy to implement, naturally incorporates the boundary conditions, and suffers no accuracy loss near boundaries or material interfaces, while post-processing with a convolution requires special kernels near the boundary [1, 26, 31, 25]. A disadvantage of the proposed method is that it requires a global operator, while the convolution with a Kernel is typically a local operator. This makes the new method unsuitable for problems that require an accurate reconstruction of the wave field throughout the entire time interval.

The paper is constructed as follows: first, we present the scalar wave equation and the corresponding classical finite element method in Section 2. In Section 3, we then introduce the new pre- and post-processing algorithm. Super-convergence rates of the classical finite element method in negative-order Sobolev norms and of the post-processed solution in the energy- and $L^{2}$-norm are derived in Section 4. In Section 5, we then demonstrate this super-convergence of the post-processed solution for numerous test cases, including cases with heterogeneous materials, unstructured meshes, and curved boundaries. Finally, we summarise our main conclusions in Section 6.

2. Wave equation and finite element discretisation

Let $\Omega\subset\mathbb{R}^{d}$ be an open bounded $d$-dimensional domain with Lipschitz boundary $\partial\Omega$. We consider the scalar wave equation given by

[TABLE]

where $(0,T)$ is the time domain, $u:\Omega\times(0,T)\rightarrow\mathbb{R}$ is the wave field, $\nabla$ is the gradient operator, $u_{0},v_{0}:\Omega\rightarrow\mathbb{R}$ are the initial wave and velocity field, $f:\Omega\times(0,T)\rightarrow\mathbb{R}$ is the source term, and $\rho,c:\Omega\rightarrow\mathbb{R}^{+}$ are strictly positive scalar fields that satisfy

[TABLE]

for some positive constants $\rho_{0}$, $\rho_{1}$, $c_{0}$, and $c_{1}$.

For the weak formulation, let $u_{0}\in H_{0}^{1}(\Omega)$, $v_{0}\in L^{2}(\Omega)$, and $f\in L^{2}(0,T;L^{2}(\Omega))$, where $L^{2}(\Omega)$ denotes the standard Lebesque space of square integrable functions on $\Omega$, $H^{1}_{0}(\Omega)$ denotes the standard Sobolev space of functions on $\Omega$ that have a zero trace on the boundary $\partial\Omega$ and have square-integrable weak derivatives, and $L^{2}(0,T;U)$, with $U$ a Banach space, is the standard Bochner space of functions $f:(0,T)\rightarrow U$, such that $\|f\|_{U}$ is square integrable on $(0,T)$. The weak formulation of (1) can then be formulated as finding u\in L^{2}\big{(}0,T;H^{1}_{0}(\Omega)\big{)}, with \partial_{t}u\in L^{2}\big{(}0,T;L^{2}(\Omega)\big{)}, \partial_{t}(\rho\partial_{t}u)\in L^{2}\big{(}0,T;H^{-1}(\Omega)\big{)}, such that $u|_{t=0}=u_{0}$, $\partial_{t}u|_{t=0}=v_{0}$, and

[TABLE]

where $H^{-1}(\Omega)$ denotes the dual space of $H^{1}_{0}(\Omega)$, $\langle\cdot,\cdot\rangle$ denotes the pairing between $H^{-1}(\Omega)$ and $H^{1}_{0}(\Omega)$, $(\cdot,\cdot)_{\rho}$ denotes the weighted $L^{2}(\Omega)$-product

[TABLE]

and $a(\cdot,\cdot)$ denotes the elliptic operator

[TABLE]

It can be shown, in a way analogous to [18, Chapter 3, Theorem 8.1], that this problem is well posed and has a unique solution.

The weak formulation can be solved with the classical finite element method. Let $\mathcal{T}_{h}$ be a tessellation of $\Omega$ into straight and/or curved simplicial elements that all have a radius smaller than $h$. The degree-$p$ finite element space $U_{h}$ can then be constructed as follows:

[TABLE]

where $\phi_{e}:\hat{e}\rightarrow e$ denotes the reference-to-physical element mapping and $\hat{U}:=\mathcal{P}_{p}(\hat{e})$ denotes the reference element space, with $\mathcal{P}_{p}(\hat{e})$ the space of polynomials on $\hat{e}$ of degree $p$ or less. The element mappings $\phi_{e}$ are usually affine mapping resulting in straight elements, but may also be higher-degree polynomial mappings resulting in curved elements that can better approximate domains with curved boundaries. The finite element method approximates the wave field by a discrete wave field $u_{h}:[0,T]:\rightarrow U_{h}$ that satisfies the initial conditions $u_{h}|_{t=0}=u_{0,h}$ and $\partial_{t}u_{h}|_{t=0}=v_{0,h}$ and that satisfies

[TABLE]

The initial values $u_{0,h},v_{0,h}\in U_{h}$ and the discrete source term f_{h}\in L^{2}\big{(}0,T;L^{2}(\Omega)\cap U_{h}\big{)} are projections of $u_{0}$, $v_{0}$, and $f$, onto the discrete space. We define $f_{h}:=\Pi_{h}f$, where $\Pi_{h}$ is the weighted $L^{2}$-projection operator given by

[TABLE]

Typically, a (weighted) $L^{2}$-projection is also used to discretise the initial solutions, although other types of projections are introduced in the next section in order to obtain higher convergence rates.

3. Pre- and post-processing

To present the pre- and post-processing algorithm, we first define the differential operator $\mathcal{L}$ by

[TABLE]

where $\mathcal{L}u\in L^{2}(\Omega)$ is well-defined when $u\in H^{1}(\Omega)$ and $\nabla\cdot c\nabla u\in L^{2}(\Omega)$. The inverse operator $\mathcal{L}^{-1}:L^{2}(\Omega)\rightarrow H_{0}^{1}(\Omega)$ is defined such that

[TABLE]

Note that $u=\mathcal{L}^{-1}f$ is the weak solution of the elliptic problem

[TABLE]

The existence of a unique weak solution follows from the fact that $a$ is coercive, which in turn follows from the Poincaré inequality and the boundedness of $c$. Also note that, for any $f\in L^{2}(\Omega)$, we have $\mathcal{L}\mathcal{L}^{-1}f=f$ and, for any $u\in H_{0}^{1}(\Omega)$ with $\mathcal{L}u\in L^{2}(\Omega)$, we have $\mathcal{L}^{-1}\mathcal{L}u=u$. We will also often use the notation $\mathcal{L}^{\alpha}$ and $\mathcal{L}^{-\alpha}$ to denote the operator $\mathcal{L}$ and $\mathcal{L}^{-1}$ applied $\alpha$ times, respectively, , e.g. $\mathcal{L}^{2}u=\mathcal{L}(\mathcal{L}u)$, $\mathcal{L}^{3}=\mathcal{L}(\mathcal{L}(\mathcal{L}u))$, etc.

The discrete version of $\mathcal{L}$, denoted by $\mathcal{L}_{h}:U_{h}\rightarrow U_{h}$, is defined such that

[TABLE]

and the inverse operator $\mathcal{L}_{h}^{-1}:U_{h}\rightarrow U_{h}$ such that

[TABLE]

Note that $u_{h}=\mathcal{L}_{h}^{-1}f_{h}$ is the standard finite element approximation of the solution $u$ to (5).

It can be shown that, if the weak solution $u$ of (1) satisfies \partial_{t}^{2}u\in L^{2}\big{(}0,T;L^{2}(\Omega)\big{)} and if f\in L^{2}\big{(}0,T;L^{2}(\Omega)\big{)}, then $u$ satisfies (1) almost everywhere and we can write

[TABLE]

while for the finite element approximation, we can write

[TABLE]

When $u$ and $f$ are sufficiently regular, we can differentiate (6) and (7) with respect to time and obtain equalities of the form

[TABLE]

for $k\geq 0$. By reordering the terms, we can obtain equalities of the form

[TABLE]

and

[TABLE]

for $k\geq 0$, where $r_{k}$, $r_{k}^{-1}$, $r_{k,h}$, and $r_{k,h}^{-1}$ are mappings defined as follows:

[TABLE]

for $k\geq 0$. By applying (8) and (9) multiple times, we can obtain the following relations:

[TABLE]

and

[TABLE]

for $k\geq 0,\alpha\geq 0$, with $r_{k}^{(0)}$ and $r_{k,h}^{(0)}$ the identity operators. Note that

[TABLE]

To improve the convergence rate of the finite element approximation by $q\geq 1$ orders, assume $u$ and $f$ are sufficiently regular and discretise the initial values as follows:

[TABLE]

To obtain an improved approximation $u_{h*}$ and $v_{h*}$ of the wave field $u$ and velocity $\partial_{t}u$ at time $T$, respectively, we compute

[TABLE]

Values of $\alpha$ and $\beta$ for $q\leq 10$ are listed in Table 1. In case of sufficient regularity, we can then obtain error estimates of the form

[TABLE]

(see Corollary 4.5), where $p$ denotes the degree of the polynomial approximation and $C$ is some positive constant that does not depend on the mesh resolution $h$. This implies, for example, that a convergence rate of order $2p$ in the energy norm can be obtained when choosing $q=p$.

From (10) and (11), it follows that the choice for the initial conditions is equivalent to setting

[TABLE]

It also follows from (10) and (11), that the post-processed solution at time $T$ is equal to the exact solution when

[TABLE]

In case of no source term, so $f=0$, the pre- and post-processed approximations simplify to

[TABLE]

To discretise the initial values $u_{0}$ and $v_{0}$, we need to respectively apply $\alpha$ and $\beta$ times the operator $\mathcal{L}_{h}^{-1}$, which means we need to respectively solve $\alpha$ and $\beta$ elliptic problems of the form in (5) numerically using the finite element method. The mesh and approximation space used to solve these elliptic problems are the same as those used for the time-stepping. To obtain the improved wave field $u_{h*}$ or velocity field $v_{h*}$ at the final time, we need to respectively apply $\alpha$ or $\beta$ times the operator $\mathcal{L}^{-1}$, which requires solving elliptic problems of the form in (5) $\alpha$ or $\beta$ times exactly. Obtaining the exact solution, however, is usually not possible. Instead, we can approximate these elliptic problems with a degree-$(p+q)$ finite element method in order to maintain the improved convergence rates.

In the end, both pre- and post-processing therefore require solving a set of at most $\alpha+\beta=q$ linear systems. These linear systems can be solved by a direct method or with an iterative method. When using an iterative method, good first approximations can be obtained from the unprocessed initial and final value, as shown in Section 5.

The pre- and post-processing algorithms presented here apply to wave problems with zero Dirichlet boundary conditions, but can also be applied to problems with Neumann boundary conditions or periodic domains when we replace $H_{0}^{1}(\Omega)$ by $H^{1}(\Omega)\cap\{1\}^{\perp}$, where $\{1\}^{\perp}:=\{u\in L^{2}(\Omega)\;|\;(u,1)_{\rho}=0\}$. Furthermore, the algorithms can also be applied to problems with non-zero boundary conditions, since the non-zero part can be absorbed in the source term $f$ and the initial conditions. In particular, let $\tilde{u}$ be any smooth function that satisfies the non-zero boundary condition and set $\tilde{f}:=\partial_{t}^{2}\tilde{u}-\rho^{-1}\nabla\cdot c\nabla\tilde{u}$. We can then write $u=\hat{u}+\tilde{u}$, where $\hat{u}$ is the solution of the wave equation (1) with zero boundary conditions, source term $\hat{f}=f-\tilde{f}$ and initial conditions $\hat{u}_{0}=u_{0}-\tilde{u}|_{t=0}$ and $\hat{v}_{0}=v_{0}-\partial_{t}\tilde{u}|_{t=0}$.

4. Error analysis

In this section, we analyse the accuracy of the classical finite element method in negative-order Solobev norms and the accuracy of the pre- and post-processing algorithms in the $L^{2}$- and energy norm. We only consider the semi-discrete case, so not the time discretisation. The outline of the error analysis is as follows: we first obtain error estimates for the classical finite element approximation of the form

[TABLE]

for $m\geq 0$ (see Theorem 4.2), where $\|\cdot\|_{0*}:=\|\cdot\|_{0}$ and $\|\cdot\|_{1*}:=\|\cdot\|_{1}$ denote the standard $L^{2}$- and $H^{1}$-norm, respectively, and $\|\cdot\|_{-m*}$, for $m>0$, denotes the adapted negative-order Sobolev-norm defined in (16). The theory is similar to that of [11], except that in [11], they require a very complex and non-standard discretisation of the initial values, while we prove in Theorem 4.2 that the estimates also hold when the initial values are discretised using a standard weighted $L^{2}$-projection.

We then show (see proof Theorem 4.4) that $\partial_{t}^{q}u$ and $\partial_{t}^{q}u_{h}$ can also be written as a solution and finite element approximation of the wave equation, respectively, and that by discretising $u_{0,h}$ and $v_{0,h}$ as in (12), the error is bounded by

[TABLE]

for $m\geq 0$.

Finally, we show that

[TABLE]

for $m\geq 0$ (see Theorem 4.4), where $u_{h*}$ and $v_{h*}$ are post-processed solutions, defined as in (13). By choosing $m=1$ or $m=0$, we then obtain bounds of the form

[TABLE]

The analysis presented in this section relies on certain regularity assumptions that typically only hold when the boundary of the domain is sufficiently smooth. We therefore do not restrict our analysis to polyhedral domains, but assume that $\Omega$ has a piecewise smooth boundary. An accurate finite element approximation will then require the use of curved elements and the leading constants in the error estimates will then also depend on the shape-regularity of these elements. We therefore define the shape-regularity of the mesh, denoted by $\gamma$, as follows:

[TABLE]

where $h_{e}$ is the diameter of $e$, $\|\cdot\|$ denotes the Euclidean norm, $D^{\alpha}:=\partial_{1}^{\alpha_{1}}\partial_{2}^{\alpha_{2}}\dots\partial_{d}^{\alpha_{d}}$ is a partial differential operator, $\partial_{i}$ is the partial derivative with respect to Cartesian coordinate $x_{i}$, and $|\boldsymbol{\alpha}|:=\alpha_{1}+\alpha_{2}+\dots+\alpha_{d}$ is the order of the operator.

In case of a straight element, $\gamma_{e}$ is proportional to the commonly used ratio $h_{e}/\rho_{e}$, where $\rho_{e}$ denotes the diameter of the inscribed sphere of $e$. For curved elements, the regularity constant $\gamma$ remains uniformly bounded for a sequence of meshes when using a suitable parametrisation, such as the one described in [16]. We present a simpler alternative to this algorithm in Section 5.5.

To simplify the error analysis, we will always assume that the mesh covers the domain exactly, so $\bigcup_{e\in\mathcal{T}_{h}}\overline{e}=\overline{\Omega}$. Furthermore, for the readability of the analysis, we will always let $C$ denote a constant that may depend on the domain $\Omega$, the regularity of the mesh $\gamma$, the spatial parameters $\rho$ and $c$, and the polynomial degree $p$, but does not depend on the mesh resolution $h$, time interval $(0,T)$, or the functions that appear in the inequality. We will also always assume that $h\leq h_{0}$, for some sufficiently small $h_{0}>0$ that does not depend on the time interval $(0,T)$ or the functions that appear in the inequality.

Let $H^{k}(\Omega)$, for $k\geq 1$, denote the standard Sobolev space of functions on $\Omega$ with square-integrable $k$th-order derivatives equipped with norm

[TABLE]

where $\|\cdot\|_{0}$ denotes the standard $L^{2}(\Omega)$-norm. Also, let $H^{k}_{0}(\Omega)$ denote the space of functions $u\in H^{k}(\Omega)$ such that the trace of the derivatives is zero on the boundary: $D^{\boldsymbol{\alpha}}u|_{\partial\Omega}=0$, for all $|\boldsymbol{\alpha}|\leq k-1$. We then define the following negative-order Sobolev norms for functions in $L^{2}(\Omega)$:

[TABLE]

We also introduce an adapted version of these negative order norms:

[TABLE]

It is proven in Section C, that, if $\rho$ and $c$ are sufficiently smooth, then

[TABLE]

Now let $U$ be a Banach space and let $L^{q}(0,T;U))$ denote the Bochner space of functions $f:(0,T)\rightarrow U$ such that $\|f(t)\|_{U}$ is in $L^{q}(0,T)$. We equip the space with norm

[TABLE]

We will use $\|\cdot\|_{\infty,k}$ as a short-hand notation for the $L^{\infty}(0,T;H^{k}(\Omega))$-norm and define $\|\cdot\|_{\infty,-k*}$, for $k\geq 1$, as follows:

[TABLE]

We will often make a regularity assumption of the following type: if $f\in H^{k}(\Omega)$, then $\mathcal{L}^{-1}f\in H^{k+2}(\Omega)$ and

[TABLE]

Such an assumption certainly holds when $\partial\Omega\in\mathcal{C}^{K+2}$, $c\in\mathcal{C}^{K+1}(\overline{\Omega})$, and $\rho\in\mathcal{C}^{K}(\overline{\Omega})$.

We can now prove super-convergence for the finite element method combined with pre- and post-processing by proving the following lemma’s and theorems.

Lemma 4.1** (error equation).**

Let $u$ be the solution of (3) and let $u_{h}$ be the solution of (4) with $f_{h}:=\Pi_{h}f$. Assume \partial_{t}^{2}u\in L^{2}\big{(}0,T;L^{2}(\Omega)\big{)}. Then the discrete error $e_{h}:=\Pi_{h}u-u_{h}$ satisfies

[TABLE]

for all $w\in U_{h}$ and almost every $t\in(0,T)$, where $\epsilon_{h}:=u-\Pi_{h}u$ is the projection error.

Proof.

By subtracting (4) from (3) and by using the assumption that \partial_{t}^{2}u\in L^{2}\big{(}0,T;L^{2}(\Omega)\big{)}, we obtain

[TABLE]

for all $w\in U_{h}$, almost every $t\in(0,T)$. Substituting $u-u_{h}=e_{h}+\epsilon_{h}$ and reordering the terms gives

[TABLE]

for all $w\in U_{h}$, almost every $t\in(0,T)$. By definition of $\Pi_{h}$, $(\partial_{t}^{2}\epsilon_{h},w)_{\rho}=0$ and $(f-f_{h},w)_{\rho}=(f-\Pi_{h}f,w)_{\rho}=0$ for all $w\in U_{h}$. The above then results in (18). ∎

Theorem 4.2**.**

Let $u$ be the weak solution of (3) and let $u_{h}$ be the degree-$p$ finite element approximation of (4), with $p\geq 1$ and with $u_{0,h}:=\Pi_{h}u_{0}$, $v_{0,h}:=\Pi_{h}v_{0}$, and $f_{h}:=\Pi_{h}f$. Assume regularity condition (17) holds for some $K\geq 0$ and assume u,\partial_{t}u\in L^{\infty}\big{(}0,T;H^{k}(\Omega)\big{)} for some $k\geq 1$ and \partial_{t}^{2}u\in L^{\infty}\big{(}0,T;L^{2}(\Omega)\big{)}. Then

[TABLE]

for any $m\geq 1$.

Proof.

We first consider the case $m=2\alpha$, with $\alpha\geq 0$. Define discrete error $e_{h}:=\Pi_{h}u-u_{h}$ and projection error $\epsilon_{h}:=u-\Pi_{h}u$ and use Lemma 4.1 to obtain

[TABLE]

for all $w\in U_{h}$, almost every $t\in(0,T)$. Set $w=\mathcal{L}_{h}^{-2\alpha}\partial_{t}e_{h}$ and use Lemma A.1 to obtain

[TABLE]

where

[TABLE]

Integrating over $(0,T^{\prime})$, with $T^{\prime}\in(0,T)$, results in

[TABLE]

From the boundedness of $\rho$ and the coercivity of $a$, it follows that

[TABLE]

To bound $R_{h}$, we derive the following:

[TABLE]

for $i=0,1$ and almost every $t\in(0,T)$. Here, the first line follows from the triangle inequality, the second line follows from Lemma A.1, the third line follows from the Cauchy–Scwarz inequality and Lemma B.7, and the last line follows from Lemma B.1 and Lemma B.3. Using this inequality, we can obtain

[TABLE]

for almost every $T^{\prime}\in(0,T)$. Since $u_{h,0}=\Pi_{h}u_{0}$ and $v_{h,0}=\Pi_{h}v_{0}$, we have $e_{h}|_{t=0}=0$ and $\partial_{t}e_{h}|_{t=0}=0$ and therefore

[TABLE]

Taking the essential supremum of (21) over all $T^{\prime}\in(0,T)$ and using (22), (23), and (24) results in

[TABLE]

Using this result, we can derive

[TABLE]

Here, we used the triangle inequality in the second line, Lemma B.7 in the third line, and (25) in the last line. Since $u-u_{h}=e_{h}+\epsilon_{h}$, (19) then follows from the above and Lemma B.3.

In case $m=2\alpha+1$, we choose $w$ in (20) as $w=\mathcal{L}_{h}^{-(2\alpha+1)}$. Using Lemma A.1, we can then derive

[TABLE]

where

[TABLE]

The rest of the proof is analogous to the proof for the case $m=2\alpha$. ∎

Lemma 4.3**.**

Let w\in L^{\infty}\big{(}0,T;L^{2}(\Omega)\big{)} and w_{h}\in L^{\infty}\big{(}0,T;U_{h}\big{)} be two arbitrary functions. Also, let $k\geq 0$ and $\alpha\geq 1$ and assume that \partial_{t}^{i}f\in L^{\infty}\big{(}0,T;L^{2}(\Omega)\big{)} for all $i\leq k+2\alpha-2$. Then

[TABLE]

Proof.

The results follow readily from induction on $\alpha$. ∎

Theorem 4.4**.**

Let $u$ be the weak solution of (3) and let $u_{h}$ be the degree-$p$ finite element approximation of (4), with $p\geq 1$ and with $f_{h}:=\Pi_{h}f$. Assume regularity condition (17) holds for some $K\geq 0$. Also, let $q\geq 0$ and $k\geq 1$ and assume that \partial_{t}^{i}u\in L^{\infty}\big{(}0,T;H^{k}(\Omega)\big{)} for all $i\leq q+1$, \partial_{t}^{q+2}u\in L^{\infty}\big{(}0,T;L^{2}(\Omega)\big{)} and \partial_{t}^{i}f\in L^{\infty}\big{(}0,T;L^{2}(\Omega)\big{)} for all $i\leq q$. If the initial conditions are discretized by

[TABLE]

and if the post-processed solution is computed by

[TABLE]

then

[TABLE]

for all $m\geq 0$.

Proof.

From the regularity of the time derivatives of $u$ and $f$, it follows that

[TABLE]

for a.e. $(x,t)\in\Omega\times(0,T)$. From the definitions of $r_{i}^{(\kappa)}$ and $r_{i,h}^{(\kappa)}$, it also follows that

[TABLE]

for $i+2\kappa\leq q+1$.

Now, define $u^{(q)}:=\partial_{t}^{q}u$ and $u^{(q)}_{h}:=\partial_{t}^{q}u_{h}$. From (28), it follows that $u^{(q)}$ is the solution of the weak formulation given in (3) when we replace $u_{0}$, $v_{0}$, and $f$, by $u^{(q)}_{0}:=\partial_{t}^{q}u|_{t=0}$, $u^{(q+1)}_{0}:=\partial_{t}^{q+1}u|_{t=0}$, and $f^{(q)}:=\partial_{t}^{q}f$, respectively, and $u^{(q)}_{h}$ is the finite element approximation given in (4) when we replace $u_{0,h}$, $v_{0,h}$, and $f_{h}$, by $u^{(q)}_{0,h}:=\partial_{t}^{q}u_{h}|_{t=0}$, $u^{(q+1)}_{0,h}:=\partial_{t}^{q+1}u_{h}|_{t=0}$, and $f^{(q)}_{h}:=\partial_{t}^{q}f_{h}=\Pi_{h}f^{(q)}$, respectively.

We now prove that, from the discretisation of $u_{h,0}$ and $v_{h,0}$, it follows that

[TABLE]

To prove this, we first consider the case $q=2\alpha$. Then $\beta=\alpha$ and we can derive

[TABLE]

Here, the first equality follows from (29b), the second equality from the discretisation of $u_{h,0}$, and the last equality from (29a). In an analogous way, we can show that

[TABLE]

The proof of (30) for the case $q=2\alpha-1$ is analogous to that for the case $q=2\alpha$. From Theorem 4.2, it then follows that

[TABLE]

for all $m^{\prime}\geq 0$.

Next, we prove that

[TABLE]

for $m\geq 0$, a.e. $t\in(0,T)$. We first consider the case $q=2\alpha$. Then $\beta=\alpha$ and we can derive

[TABLE]

Here, the first equality follows from the definition of $u_{h*}$, the second equality follows from (29), the third equality follows from Lemma 4.3, and the last equality follows from the definitions of the adapted negative-order Sobolev norm. In an analogous way, we can derive

[TABLE]

Together, these last two equalities result in (32) for the case $q=2\alpha$. The proof for the case $q=2\alpha-1$ is analogous to that for $q=\alpha$. Inequality (27) then follows immediately from (31) and (32) when we set $m^{\prime}=q+m$. ∎

By taking $m=0$ and $m=1$ in Theorem 4.4, we immediately obtain the following.

Corollary 4.5**.**

Let $u$ be the weak solution of (3), let $u_{h}$ be the degree-$p$ finite element approximation of (4), with $p\geq 1$ and with $f_{h}$, $u_{0,h}$, $v_{0,h}$, $u_{h*}$, and $v_{h*}$ chosen as in Theorem 4.4. If the assumptions in Theorem4.4 hold for a certain $q\geq 0$ and for $k=p+1$ and $K=p-1$, then

[TABLE]

and

[TABLE]

Remark 4.6**.**

While the leading constant in Theorem 4.4 depends on the shape-regularity of the elements, the theory is not restricted to (quasi-)uniform meshes but holds for general unstructured meshes.

5. Implementation and numerical tests

5.1. Time discretisation

In the error analysis, we considered exact integration in space and time. In practice, however, we also need to discretise in time and use numerical integration to evaluate the spatial integrals and lump the mass matrix.

To discretise in time, we first rewrite the finite element formulation given in (4) as a system of ODE’s. To do this, we use nodal basis functions. Let $\mathcal{Q}_{h}=\{\mathbf{x}_{i}\}_{i=1}^{N}$ be a set of nodes of the form

[TABLE]

where $\hat{\mathcal{Q}}_{I}$ are the nodes on the reference element, and let $\{w_{i}\}_{i=1}^{N}$ be the nodal basis functions that span the discrete space $U_{h}$ and satisfy $w_{i}(\mathbf{x}_{j})=\delta_{ij}$, for $i,j=1,\dots,N$, with $\delta$ the Kronecker delta. Also, for any $f\in\mathcal{C}^{0}(\overline{\Omega})$, define the vector $\underline{f}\in\mathbb{R}^{N}$ as $\underline{f}_{i}:=f(\mathbf{x}_{i})$ for $i=1,\dots,N$. The finite element formulation can then be written as finding $\underline{u_{h}}:[0,T]\rightarrow\mathbb{R}^{n}$, such that $\underline{u_{h}}|_{t=0}=\underline{u_{0,h}}$, $\partial_{t}\underline{u_{h}}|_{t=0}=\underline{v_{0,h}}$, and

[TABLE]

where $L_{h}:=M^{-1}A$ and $M,A\in\mathbb{R}^{N\times N}$ are the mass matrix and stiffness matrix, respectively, given by

[TABLE]

For the time discretisation, let $\Delta t$ denote the time step size, let $t_{n}:=n\Delta t$, and let $\underline{u}_{h}^{n}$ denote the approximation of $\underline{u_{h}}(t_{n})$. In order to maintain an order-$2p$ convergence rate, we use an order-$2p$ Dablain scheme [10] for time-stepping, which is given by

[TABLE]

where $\underline{\partial_{t}^{k}u}_{h}^{n}$ is recursively defined by

[TABLE]

so the approximation $\underline{u}_{h}^{n+1}$ is computed using the two previous approximations $\underline{u}_{h}^{n-1}$ and $\underline{u}_{h}^{n}$. For the first time step, the computations are as follows:

[TABLE]

with $\underline{u}_{h}^{0}:=\underline{u_{0,h}}$ and $\underline{\partial_{t}u}_{h}^{0}:=\underline{v_{0,h}}$. This scheme is commonly used for wave propagation modelling and has the advanatge that it only requires $p$ stages to obtain an order-$2p$ convergence rate.

Dablain’s scheme only computes the displacement field. A second-order approximation of the velocity $\underline{\partial_{t}u_{h}}$ at some time $t_{n}$ can be obtained by

[TABLE]

In case $\underline{u}_{h}^{n+1}=\underline{u_{h}}(t_{n+1})$, $\underline{u}_{h}^{n-1}=\underline{u_{h}}(t_{n-1})$, and $\underline{\partial_{t}u}_{h}^{n}=\underline{\partial_{t}u_{h}}(t_{n})$, the Taylor approximation of $\underline{v}_{h}^{n,2}$ around $t_{n}$ is given by

[TABLE]

Using this expansion, we can obtain higher-order approximations of the velocity. For example, fourth- and sixth-order approximations are given by

[TABLE]

Since we are only interested in the order-$2p$ accurate velocity field at the final time slot $n=N_{T}:=T/\Delta t$, we only need to do this computation once.

For stability, the time step size should be sufficiently small:

[TABLE]

where $c_{p}=4,12,7.57$ [12], for $p=1,2,3$, respectively, and $\sigma_{\max}(L_{h})$ denotes the largest eigenvalue of $L_{h}$. The largest eigenvalue can be bounded by the largest eigenvalues of the element matrices [14]:

[TABLE]

where $M_{e},A_{e}$ denote the mass- and stiffness matrix of element $e$, respectively. In the numerical tests, we will always set $\Delta t=0.9\Delta t_{\max}$, with $\sigma_{\max}(L_{h})$ estimated using the above.

When pre-processing, the initial discrete values $\underline{u_{0,h}}=\underline{u}_{h}^{0}$ and $\underline{v_{0,h}}=\underline{\partial_{t}u}_{h}^{0}$ are computed by

[TABLE]

and by recursively solving

[TABLE]

The initial values $\partial_{t}^{q}u|_{t=0}$ and $\partial_{t}^{q+1}u|_{t=0}$ can be obtained from $u_{0}$ and $v_{0}$ by computing (8a) for $k=0,1,\dots,q-1$.

When post-processing, we need to apply operators of the form $r_{k}^{(-1)}$ and therefore apply the operator $\mathcal{L}^{-1}f$, which maps a function $f$ to the exact solution of the elliptic problem given in (5). Since it is usually not possible to solve the elliptic problem exactly, we approximate $\mathcal{L}^{-1}$ by $\mathcal{L}_{h*}^{-1}$, which maps a function $f$ to a finite element approximation of the elliptic problem. For this finite element approximation, we use the same mesh as for the time-stepping but with a degree-$(p+q)$ finite element space.

Let $U_{h*}$ denote the higher-degree finite element space, $\mathcal{Q}_{h*}=\{\mathbf{x}_{i*}\}_{i=1}^{N_{*}}$ the corresponding nodes, $M_{*}$, $A_{*}$, and $L_{h*}:=M^{-1}_{*}A_{*}$ the corresponding matrices, and $\underline{f}_{*}$ the vector of values of a function $f$ at the nodes $\mathbf{x}_{i*}$. Also, let $P\in\mathbb{R}^{N_{*}\times N}$, defined by $P_{ij}:=w_{j}(\mathbf{x}_{i*})$, be the matrix that maps the degrees of freedom of the degree-$p$ finite element space to the degrees of freedom of the degree-$(p+q)$ finite element space. The post-processed wave field $\underline{u}_{h*}^{n}$ and velocity field $\underline{v}_{h*}^{n}:=\underline{\partial_{t}u}_{h*}^{n}$ can then be computed by

[TABLE]

and by recursive solving

[TABLE]

The discrete time derivatives $\underline{\partial_{t}^{q}u}_{h}^{n}$ and $\underline{\partial_{t}^{q+1}u}_{h}^{n}$ can be obtained from $\underline{u}_{h}^{n}$ and $\underline{\partial_{t}u}_{h}^{n}:=\underline{v}_{h}^{n,2p}$ by recursively computing (36) for $k=0,1,\dots,q-1$. Again, we only need to do these post-processing steps at the final time slot $n=N_{T}$.

5.2. Quadrature rules and mass lumping

To compute the spatial integrals, we use a quadrature rule that consists of a set of points $\hat{\mathcal{Q}}_{M}$ on the reference element and a set of weights $\{\omega_{\hat{\mathbf{x}}}\}_{\hat{\mathbf{x}}\in\hat{\mathcal{Q}}_{M}}$. The integral over the reference element is then approximated as follows:

[TABLE]

We can write integrals over the physical element as integrals over the reference element using the relation

[TABLE]

where $\mathbf{J}_{e}:=\hat{\nabla}\phi_{e}$ denotes the Jacobian of the element mapping, $\hat{\nabla}$ is the gradient operator of the reference space, and $|\mathbf{J}_{e}|:=|\mathrm{det}(\mathbf{J}_{e})|$.

A major drawback of the classical finite element method for the wave equation is that the mass matrix is not strictly diagonal. Since time stepping requires computing terms of the form $L_{h}\underline{u}=M^{-1}A\underline{u}$ at each time step, a non-diagonal mass matrix would require solving a large system of equations at each time step. In practice, the mass matrix is therefore lumped into a diagonal matrix, so that the system of equations becomes trivial to solve. A lumped mass matrix can be obtained by placing the nodes of the basis functions $\hat{\mathcal{Q}}_{I}$ at the quadrature points for the mass matrix, so by setting $\hat{\mathcal{Q}}_{I}=\hat{\mathcal{Q}}_{M}$. To obtain 1D mass-lumped elements, we can use Gauss-Lobatto points. This can be extended to quadrilateral and hexahedral mass-lumped elements by using tensor-product basis functions. The resulting scheme is known as the spectral element method [24, 28, 15]. For linear triangular and tetrahedral elements, we can place the nodes at the vertices and use a Newton–Cotes integration rule. For higher-order triangular and tetrahedral elements, however, we need to enrich the element space with higher-degree bubble functions and use special quadrature rules to maintain stability and accuracy after mass-lumping [7, 6, 22, 4, 23, 19, 8, 13].

In this paper, we test using the 1D spectral elements, the linear mass-lumped triangular element, and the quadratic and cubic mass-lumped triangular element presented in [7]. The element space of the quadratic mass-lumped triangular element is given by

[TABLE]

where $b:=\hat{x}_{1}\hat{x}_{2}\hat{x}_{3}$ is the degree-3 bubble function, with $\hat{x}_{i}$ the barycentric coordinates. The nodes of this element are placed at the 3 vertices, the midpoint of the 3 edges, and the centre of the triangle. The space of the cubic mass-lumped triangular element is given by

[TABLE]

The nodes of this element are placed at the three vertices, the 6 points on the edges with barycentric coordinates $\alpha$, $1-\alpha$, and [math], and at the three interior points with barycentric coordinates $\beta$, $\beta$, and $1-2\beta$, where

[TABLE]

An illustration of these elements is given in Figure 1.

5.3. Test with sharp contrasts in material parameters

First, we test the finite element method with and without pre- and post-processing on a 1D periodic domain with sharp contrasts in material parameters and using a mesh with sharp contrasts in the element size. The test is similar to those in [27], except that, unlike in [27], super-convergence is observed on the entire domain, also near the material interfaces and sharp contrasts in the element size.

Let $\Omega:=(0,5)$ be the periodic domain and let the parameters $\rho$ and $c$ be given by

[TABLE]

The exact solution is given by the travelling wave $u(x,t):=\sin(2\pi(X(x)-t))$, where

[TABLE]

The simulated time interval is $(0,T)$ with $T=10$ and the initial conditions are obtained from the exact solution.

We solve the wave equation numerically using degree-$p$ spectral elements for the time-stepping scheme and degree-$2p$ spectral elements for the order-$q$ improving pre- and post-processing steps. In case $q=0$, no pre- and post-processing is applied. We use $N$ elements per wavelength, so $\Delta x=\frac{1}{N}$ in $(0,1)$ and $\Delta x=\frac{4}{N}$ in $(1,5)$ and therefore we have a sharp contrast in mesh size at $x=0$ and $x=1$. An illustration of the mesh is given in Figure 2

The relative errors in the weighted energy norm and weighted $L^{2}$-norm are defined as follows:

[TABLE]

Results for different elements and different pre/post-processing schemes are given in Tables 2 and 3. The errors are computed using the quadrature rule of the degree-$2p$ spectral element method. The results confirm that the convergence rate is of order $q+p$ in the energy norm and $q+p+1$ in the $L^{2}$-norm. In case $p=3$ and $q=1$, this scheme even seems to converge with order 6 in the $L^{2}$-norm. However, this higher convergence rate only appears in the 1D case and not in the 2D case as we will show in the next subsection.

Figure 3 also shows the error of the unprocessed and post-processed finite element approximation for $N=10$ and $p=3$. The figure illustrates that the error of the unprocessed approximation is highly oscillatory, while the error of the post-processed solution is much smoother.

Since there are strong contrasts in the material parameters, it is not immediately obvious that the regularity assumption given in (17) holds. However, by using the piecewise-linear coordinate transformation $x\rightarrow X$, given in (39), the problem can be rewritten as a wave propagation problem for a periodic domain with constant material parameters.

5.4. Test using unstructured triangular meshes

Next, we test the pre- and post-processing methods on a squared domain with non-constant material parameters and zero Dirichlet boundary conditions using unstructured triangular meshes. The test is similar to those in [21], except that, unlike in [21], super-convergence is still clearly observed when using unstructured meshes.

Let $\Omega:=[0,1]^{2}$ be the spatial domain and $(0,T)$ the time domain with final time $T=50$. Introduce distortions of the Cartesian coordinates, given by

[TABLE]

and define the spatial parameters

[TABLE]

where $\dot{X}$ and $\dot{Y}$ denote the derivatives of $X$ and $Y$, respectively. The wave propagation speed $\sqrt{c/\rho}$ then takes values between $(1+\frac{1}{10}\pi)^{-1}\approx 0.76$ and $(1-\frac{1}{10}\pi)^{-1}\approx 1.46$. The exact solution is given by

[TABLE]

with source term

[TABLE]

where $\phi:=0.5$ is an arbitrary phase shift added to avoid the initial displacement or velocity field to be completely zero. The initial conditions are obtained from the exact solution.

We consider 5 different meshes where each subsequent mesh has a resolution twice as fine as the previous mesh. An illustration of the first three meshes is given in Figure 4. We test the standard linear mass-lumped triangular element and the quadratic and cubic mass-lumped triangular elements presented in [7]. For post-processing, we use standard degree-$2p$ Lagrangian elements combined with a degree-$4p$ accurate quadrature rule taken from [32] for evaluating the integrals. This last quadrature rule is also used to evaluate the errors.

Results for different elements and different pre/post-processing schemes are given in Tables 4 and 5. In all cases, the convergence rate is of order $q+p$ in the energy norm and $q+p+1$ in the $L^{2}$-norm. Figure 5 also shows the error of the unprocessed and post-processed finite element approximation for mesh 1 and $p=3$. The figure illustrates again that the error of the unprocessed approximation is highly oscillatory, while the error of the post-processed solution is relatively smooth.

Since the boundary of the squared domain is not smooth, it is again not immediately obvious that the regularity assumption given in (17) holds. However, by using the smooth coordinate transformations $(x,y)\rightarrow(X,Y)$, given in (40), the problem can be rewritten as a wave propagation problem on a unit square with constant material parameters. The regularity assumption (17) then follows from Fourier analysis.

5.5. Test on a domain with a curved boundary

Finally, we test the pre- and post-processing algorithms on a circular domain with zero Dirichlet boundary conditions. Let $\Omega$ be the unit circle , let $(0,T)$ be the time domain with $T=50$, and let $\rho=c=1$ be constant. The exact solution is given by

[TABLE]

where $r\in[0,1)$ and $\theta\in[0,2\pi)$ are the polar coordinates, $\phi:=0.5$ is an arbitrary phase shift, $J_{2}$ is the Bessel function of order 2, and $\kappa\approx 5.135622301840683$ is the first positive root of $J_{2}$. The initial conditions are obtained from the exact solution.

We test again on 5 different meshes where each subsequent mesh has a resolution twice as fine as the previous mesh. An illustration of the first three meshes is given in Figure 6. To obtain a higher accuracy, we replace the straight elements at the boundary by curved elements. To parametrise these curved elements, we use degree-$2p$ hyperparametric element mappings. First, we place $2p+1$ nodes at equal intervals on every edge of a boundary element, with the outer nodes lying on the vertices. We do something similar for the edges of the reference element. We then rescale the coordinates of the nodes on boundary edges such that they lie on the boundary of the unit circle. Let $\{\hat{\mathbf{x}}_{i}\}_{i=1}^{6p}$ denote the nodes at the boundary of the reference element, and let $\{\mathbf{x}_{i}\}_{i=1}^{6p}$ denote the corresponding nodes at the boundary of a curved element. We then parameterise the curved element with the element mapping $\phi_{e}$ defined such that $\phi_{e}(\hat{\mathbf{x}}_{i})=\mathbf{x}_{i}$ for $i=1,\dots,6p$. To obtain this element mapping, we interpolate using the following set of polynomials of degree $2p$ and less:

[TABLE]

where $\hat{x}_{1},\hat{x}_{2},\hat{x}_{3}$ denote the barycentric coordinates of the reference element.

For the time-stepping and pre- and post-processing, we use the same elements and quadrature rules as in the previous subsection. The results for different elements and different pre/post-processing schemes are given in Table 6. Again, the results confirm that using pre- and post-processing can result in a convergence rate of order $2p$ instead of order $p$ in the energy norm.

5.6. Pre- and post-processing using an iterative method

In the previous subsections, we applied pre- and post-processing in combination with a direct solver. In practice, the number of degrees of freedom might become too large to efficiently use a direct solver. We therefore also test the pre- and post-processing method using an iterative method.

We consider again the test on the heterogeneous squared domain using the unstructured triangular meshes of Section 5.4. For pre- and post-processing, we need to recursively compute $\underline{\partial_{t}^{k}u}_{h}^{0}$ and $\underline{\partial_{t}^{k}u}_{h*}^{N_{T}}$ using (37) and (38), respectively, which requires solving linear systems of equations of the form $A\underline{x}=\underline{y}$ and $A_{*}\underline{x}_{*}=\underline{y}_{*}$. We solve these systems using the conjugate gradient method with $N_{it}$ iterations. As a preconditioner, we use a diagonal matrix computed by taking row sums of the absolute values of $A$ and $A_{*}$, and as an initial guess for $\underline{\partial_{t}^{k}u}_{h}^{0}$ and $\underline{\partial_{t}^{k}u}_{h*}^{N_{T}}$, we use $\underline{\partial_{t}^{k}u|_{t=0}}$ and $P\underline{\partial_{t}^{k}u}_{h}^{N_{T}}$, respectively.

The accuracy of the order-$p$ improving pre- and post-processing scheme using this iterative solver with different numbers of iterations $N_{it}$ is shown in Table 7, where $N_{it}=0$ means no pre- and post-processing is applied at all and $N_{it}=\infty$ means a direct solver was used. The table shows that $10$ iterations already reduces the error by an order of magnitude, $100$ iterations reduce the error by 2-3 orders of magnitude when using higher-degree elements, and $1000$ iterations reduce the error by 4-5 orders when using degree-3 elements on finer meshes.

6. Conclusion

We presented a new pre- and post-processing method to enhance the accuracy of the finite element approximation of linear wave problems. We proved that, by applying at most $q$ processing steps at only the initial and final time, the convergence rate of the finite element approximation can be improved from order $p+1$ to order $p+1+q$ in the $L^{2}$-norm, and from order $p$ to order $p+q$ in the energy norm, both up to a maximum of order $2p$. Each processing step corresponds to solving a linear system that can be solved efficiently with a direct solver or an iterative method and for the latter, good first guesses are provided by the unprocessed initial and final values. Numerical experiments showed that the pre- and post-processing steps can reduce the magnitude of the error by several orders, even when using an iterative method with only a small number of iterations. The experiments also confirmed that a convergence rate of order $2p$ in the energy norm is obtained, even when using unstructured meshes or meshes with sharp contrasts in element size and when solving wave propagation problems on domains with curved boundaries or strong heterogeneities.

Appendix A Properties of the spatial operator

Lemma A.1**.**

The operators $\mathcal{L}^{-1}$ and $\mathcal{L}_{h}^{-1}$ satisfy the following symmetry properties:

[TABLE]

Proof.

Using the definition of $\mathcal{L}^{-1}$, we can derive the first two equalities as follows:

[TABLE]

The results for $\mathcal{L}_{h}^{-1}$ can be derived in an analogous way. ∎

Appendix B Finite element approximation properties

Lemma B.1**.**

Let $p\geq 1$ denote the degree of the finite element space and let $u\in H^{k}(\Omega)$ for some $k\geq 0$. Then

[TABLE]

Proof.

In [2], it is shown that these inequalities hold if we replace $\Pi_{h}$ by their quasi-interpolation operator $\Pi_{h}^{0}$, also when using curved elements. Inequality (42a) then follows from the fact that $\Pi_{h}$ minimises the approximation error in the weighted $L^{2}$-norm and inequality (42b) then follows from the inverse inequality. ∎

Remark B.2**.**

The leading constants in Lemma B.1 depend on the mesh regularity $\gamma$ as defined in (15). For a sequence of meshes with curved elements, this regularity constant only remains bounded for an appropriate parameterisation of the curved elements, such as, for example, the parametrisation presented in [16].

Lemma B.3**.**

Assume regularity condition (17) holds for some $K\geq 0$ and let $u\in H^{k}(\Omega)$, with $k\geq 0$. Then, for all $\alpha\geq 0$, we have

[TABLE]

where $p\geq 1$ denotes the degree of the finite element space.

Proof.

If $\alpha=0$, then the result follows immediately from Lemma B.1. Now, let $\alpha\geq 1$. Inequality (43a) immediately follows from the following:

[TABLE]

where $I$ denotes the identity operator and where the first line follows from the boundedness of $\rho$, the second line from Lemma A.1, the third line from the definition of $\Pi_{h}$, the fourth line from the Cauchy–Schwarz inequality and the boundedness of $\rho$, and the last line from Lemma B.1 and the regularity assumption.

To prove (43b), we start as follows:

[TABLE]

where the first line follows from the coercivity of $a$ and the second line from Lemma A.1 and the definition of $\mathcal{L}^{-1}$. For the rest of the proof, we can proceed in a way analogous to before. ∎

Lemma B.4**.**

Assume regularity condition (17) holds for some $K\geq 0$ and let $f\in H^{k}(\Omega)$, with $k\geq 0$. Then

[TABLE]

where $p\geq 1$ denotes the degree of the finite element space.

Proof.

This is a standard result of finite element approximations for elliptic problems. Define $u:=\mathcal{L}^{-1}f$ and $u_{h}:=\mathcal{L}_{h}^{-1}\Pi_{h}f$. From the definitions of $\mathcal{L}^{-1}$, $\mathcal{L}_{h}^{-1}$, and $\Pi_{h}$, it follows that

[TABLE]

and from the regularity assumption that $\|u\|_{\min(k+2,K+2)}\leq C\|f\|_{\min(k,K)}$. The lemma then follows from the coercivity and boundedness of $a$, Cea’s Lemma, and Lemma B.1. ∎

Lemma B.5**.**

Assume regularity condition (17) holds for some $K\geq 0$ and let $f\in H^{k}(\Omega)$, with $k\geq 0$. Then

[TABLE]

for all $\alpha\geq 0$, where $p\geq 1$ denotes the degree of the finite element space.

Proof.

Define $u:=\mathcal{L}^{-1}f$ and $u_{h}:=\mathcal{L}_{h}^{-1}\Pi_{h}f$. We first prove (45a) by deriving

[TABLE]

where $I$ denotes the identity operator and where we used the coercivity of $a$ in the first line, Lemma A.1 in the second line, the definitions of $u$ and $u_{h}$ in the third line, the boundedness of $c$ and the Cauchy–Schwarz inequality in the fourth line, and Lemma B.1 and the regularity assumption in the last line.

To prove (45b), we first derive

[TABLE]

where we used the boundedness of $\rho$ in the first line and Lemma A.1 and the definition of $\mathcal{L}^{-1}$ in the second line. For the rest of the proof, we can proceed in a way analogous to before. ∎

Corollary B.6**.**

Assume regularity condition (17) holds for some $K\geq 0$ and let $f\in H^{k}(\Omega)$, with $k\geq 0$. Furthermore, let $\mathcal{S}$ be an operator of the form

[TABLE]

with $n\geq 1$, $\beta_{i}\geq 0$ for $i=1,\dots,n+1$. Define $\alpha:=n+\beta_{1}+\beta_{2}+\cdots+\beta_{n+1}$. Then

[TABLE]

where $p\geq 1$ denotes the degree of the finite element space.

Proof.

The asserted follows from Lemma B.5 and can be proven by induction on $n$. ∎

Lemma B.7**.**

Assume regularity condition (17) holds for some $K\geq 0$ and let $f\in H^{k}(\Omega)$, with $k\geq 0$. Furthermore, let $\mathcal{S}_{1},\mathcal{S}_{2}$ be two operators of the form

[TABLE]

with $n\geq 0$, $\beta_{i}\geq 0$ for $i=1,\dots,n+1$. Define $\alpha:=n+\beta_{1}+\beta_{2}+\cdots+\beta_{n+1}$. If $\alpha\geq 1$, then

[TABLE]

where $p\geq 1$ denotes the degree of the finite element space.

Proof.

Replace the terms $\mathcal{L}_{h}^{-1}\Pi_{h}$ by $(\mathcal{L}_{h}^{-1}\Pi_{h}-\mathcal{L}^{-1})+\mathcal{L}^{-1}$ in both $\mathcal{S}_{1}$ and $\mathcal{S}_{2}$. The expansions of $\mathcal{S}_{1}$ and $\mathcal{S}_{2}$ then both consist of the term $\mathcal{L}^{-\alpha}$ and several terms of the form of (46). The difference $\mathcal{S}_{1}-\mathcal{S}_{2}$ can then be written purely in terms of operators of the form of (46). The Lemma then follows from the triangle inequality and Corollary B.6. ∎

Remark B.8**.**

The operators $\mathcal{S}_{1}$ and $\mathcal{S}_{2}$ are sequences of $\alpha$ operators, where each operator is either $\mathcal{L}^{-1}$ or $\mathcal{L}_{h}^{-1}\Pi_{h}$.

Appendix C Properties of negative-order Sobolev norms

Lemma C.1**.**

Let $u\in L^{2}(\Omega)$ and assume $\rho\in\mathcal{C}^{k}$ and $c\in\mathcal{C}^{k-1}$ for some $k\geq 0$. Then

[TABLE]

Proof.

Let $w\in H_{0}^{k}(\Omega)$ and suppose $k=2\alpha$. From the regularity of $\rho$ and $c$, it follows that $\mathcal{L}^{\kappa}(\rho^{-1}w)\in H_{0}^{1}(\Omega)$, for $\kappa\leq\alpha-1$, and $\mathcal{L}^{\alpha}(\rho^{-1}w)\in L^{2}{\Omega}$ and therefore $\mathcal{L}^{-\alpha}\mathcal{L}^{\alpha}(\rho^{-1}w)=\rho^{-1}w$. We can then derive

[TABLE]

where the first equality follows from the definition of $(\cdot,\cdot)_{\rho}$, the second equality follows from $\mathcal{L}^{-\alpha}\mathcal{L}^{\alpha}(\rho^{-1}w)=\rho^{-1}w$ and Lemma A.1, and the final inequality follows from the Cauchy–Schwarz inequality and the regularity of $\rho$ and $c$. In case $k=2\alpha+1$, we can, in an analogous way, derive

[TABLE]

Inequality (49) then follows immediately from (50) and (51). ∎