Divergence-conforming discontinuous Galerkin finite elements for Stokes eigenvalue problems

In this paper, we present a divergence-conforming discontinuous Galerkin finite element method for Stokes eigenvalue problems. We prove a priori error estimates for the eigenvalue and eigenfunction errors and present a robust residual based a posteriori error estimator. The a posteriori error estimator is proven to be reliable and (locally) efficient in a mesh-dependent velocity-pressure norm. We finally present some numerical examples that verify the a priori convergence rates and the reliability and efficiency of the residual based a posteriori error estimator.


Introduction
In fluid mechanics, eigenvalue problems are of great importance because of their role for the stability analysis of fluid flow problems.Hence, the development of numerical methods for the Stokes problem, as a model for incompressible fluid flow, is of great interest.For example in [19], several stabilized finite element methods for the Stokes eigenvalue problem are considered by Huang et al.A finite element analysis of a pseudo stress formulation for the Stokes eigenvalue problem is proposed by Meddahi et al. [27].
Currently, there are only very few results on the a posteriori error analysis for the Stokes eigenvalue problem available in the literature.An a posteriori error analysis based on residual a posteriori error estimators for the finite element discretization of the Stokes eigenvalue problem is proposed by Lovadina et al. [26].Some superconvergence results and the related recovery type a posteriori error estimators for the Stokes eigenvalue problem is presented by Liu et al. [25] based on a projection method.In [2], Armentano et al. introduced a posteriori error estimators for stabilized low-order mixed finite elements and in [16], Han et al. presented a residual type a posterior error estimator for a new adaptive mixed finite element method for the Stokes eigenvalue problem.In [18], Huang presents a posteriori lower and upper eigenvalue bounds for the Stokes eigenvalue problem for two stabilized finite element methods based on the lowest equal-order finite element pair.Recently, we have developed an a posteriori error analysis for the Arnold-Winther mixed finite element method of the Stokes eigenvalue problem in [12] using the stress-velocity formulation.
Cockburn et al. [9,10] derived the divergence-conforming discontinuous Galerkin finite element method.In [17], Houston et al. presented an a posteriori error estimation for mixed discontinuous Galerkin approximations of the Stokes problem.Kanschat et al. [22] proposed a posteriori error estimates for divergence-free discontinuous Galerkin approximations of the Navier-Stokes equations.Multigrid methods for H div -conforming discontinuous Galerkin (H div -DG) finite element methods for the Stokes equations are proposed by Kanschat et al. [21].Recently, Kanschat et al. [23] presented the relation between the H div -DG finite element method for the Stokes equation and the C 0 interior penalty finite element method for the biharmonic problem.
In this paper, we introduce an H div -DG finite element method for Stokes eigenvalue problems.We derive a priori error estimates for the eigenvalue and eigenfunction errors.We present a robust a posteriori error analysis of the H div -DG finite element method and derive upper and local lower bounds for the velocity-pressure error which is measured in terms of the mesh-dependent DG norm.The proposed a posteriori error estimator is robust in the sense that the ratio of upper and lower bounds is independent of the viscosity coefficient and the local mesh size.
For simplicity of the presentation we restrict the analysis to the case of a simple eigenvalue λ.The results can be applied to multiple eigenvalues by extending the given analysis to subspaces of eigenvectors that belong to the same multiple eigenvalue.The a posteriori error estimator can be extended to multiple eigenvalues in that the squared sum over all estimators of discrete eigenfunctions approximating the same multiple eigenvalue provides an upper bound of the eigenvalue error up to higher order terms.
The paper is organised as follows: the necessary notation and the H div -DG formulation of the Stokes eigenvalue problem is presented in Section 2. In Section 3, the a priori error analysis is discussed.The a posteriori error analysis is developed in Section 4. Finally, Section 5 is devoted to present some numerical results for uniform and adaptive mesh refinement.Let H s (ω) be the standard Sobolev space with the associated norm • s,ω for s ≥ 0. In case of ω = Ω, we use • s instead of • s,Ω .Let H −s (ω) := (H s (ω)) * be the dual space of H s (ω).Now we extend the definitions to vector and matrix-valued functions.Let H s (ω) = H s (ω; R 2 ) and H s (ω, R 2×2 ) be the Sobolev spaces over the set of 2-dimensional vector and 2 × 2 matrix-valued function, respectively.The symbols and are used to denote bounds which are valid up to positive constants independent of the local mesh size.

Notation
Throughout the paper, we consider the following spaces L 2 0 (Ω), H div (Ω), H div 0 (Ω) and ) which are defined as follows:

Weak formulation of the Stokes eigenvalue problem
Let u be the velocity, p the pressure, ν > 0 the (constant) viscosity, and Ω ⊂ R 2 be a bounded, and connected Lipschitz domain.Consider the velocity-pressure formulation of the Stokes eigenvalue problem: find an eigentripel (u, p, λ), u = 0, such that with the compatibility relation The weak formulation of the Stokes eigenvalue problem (1) reads: find (u, p, λ) We can formulate the weak formulation of (2) in a global form as: find (u, p, λ) where

Meshes, trace operators and discrete spaces
We suppose that the domain Ω is decomposed by a subdivision T h into a mesh of shape-regular rectangular cells K. Let E h denote the set of edges, E i h the set of interior edges, and E ∂ h the set of boundary edges of T h .We restrict ourself to one-irregular meshes T h in which each interior edge E ∈ E i h may contain at most one hanging node in the midpoint of E. For a given mesh T h , the notions of broken spaces for the continuous and differentiable function spaces are denoted as C 0 (T h ) and H s (T h ) which are the spaces such that the restriction to each mesh cell K ∈ T h is in C 0 (K) and H s (K), respectively.
Let K ± ∈ T h be two mesh cells which share a common edge The traces of functions v ∈ C 0 (T h ) on E from K ± are defined as v ± , respectively.Then the sum operator is defined as Let n ± be the unit outward normal vector to K ± , respectively.Then the sum operator turns into the jump operator, such that for For boundary edges E = K + ∩ ∂Ω we set [[v]] = v + and with ∇ h we denote the local application of the gradient (∇ h v)| K = ∇(v| K ) on each K ∈ T h .
We define Q k (K), Q k (K) d and Q k (K) d×d as the space of scalar, vector and tensor valued polynomials on K of partial degree at most integer k ≥ 1.
Choose V h as a discrete subspace of H div 0 (Ω) as where RT k (K) := P k+1,k (K) × P k,k+1 (K) is the Raviart-Thomas space of degree k ≥ 1, where P r,s (K) denotes the space of the polynomial functions on K of degree at most r > 0 in x 1 and at most s > 0 in x 2 .Moreover, let Q h be the discrete space of L 2 0 (Ω) such that An important property of the pair V h × Q h is as follows: on the meshes considered, see [9] for more details.As a consequence we have that the discrete velocity field u h is exactly divergence free.
Remark 2.1.The inf-sup stability of discretizations with hanging nodes using Raviart-Thomas finite elements is in part still an open question.In [29], there exists a stability proof only for the pair RT k × Q k defined in ( 4) and ( 5) with k ≥ 2 for quadrilaterals with one-irregular meshes.However, we conjecture from our computational results that stability also holds for k = 1.Moreover, the stability result for the divergence-free elements proposed in [10] is not available for triangles with hanging nodes.On the other hand, locally refined triangular meshes without hanging nodes can be obtained using bisection.The results below are all to be read in view of the restrictions cited in this remark.
Remark 2.2.The analysis of this paper also applies directly to divergence-free BDM k /P k−1 finite elements on regular triangular meshes.

H div -DG formulation for the Stokes eigenvalue problem
The discrete weak formulation of problem (1) reads: find where Here, a h (•, •) is the bilinear form defined as where the interior face terms a i p (u, v), a i c (u, v) and Nitsche terms are defined as Here, h E is the length of the edge E and γ is the penalty parameter which is chosen sufficiently large to guarantee the stability of the DG formulation, see for instance [3].
Finally, we introduce the following mesh-dependent DG velocity-pressure norm where

A priori error analysis
Our main aim is to show that the approximated eigenvalues and eigenfunctions of the H div -DG finite element formulation of the Stokes eigenvalue problem converge to the solution of the corresponding spectral problem which comes to apply the classical spectral approximation theory presented in [4,28] using results of the a priori error analysis of the associated source problem that we recall here for completeness.

Numerical analysis of the source problem
This section is devoted to discuss the source problem and to recall its essential stability and convergent results.
Consider the source problem with the right hand side f ∈ L 2 (Ω) The variational formulation of the Stokes source problem reads: find Due to the continuous inf-sup condition inf the variational formulation ( 11) is well-posed [8,15].The H div -DG finite elements of the Stokes source problem reads: find From [10,22,23], we have that the bilinear form a h (•, •) is bounded and elliptic uniformly in h on V h equipped with the norm |||•|||.Furthermore, the velocity-pressure pair V h × Q h is inf-sup stable and satisfies inf for a constant β independent of h.Hence, the weak formulation ( 13) has a unique discrete solution, which admits the following stability estimate and due to ∇ • V h ⊂ Q h the discrete velocity u f h is exactly divergence-free.From [10,14,22], we have the following a priori estimates.Theorem 3.1.Let u ∈ H s (Ω) and p ∈ H s−1 (Ω) for some s ∈]3/2, k + 1] be solutions of the continuous problem (11) that satisfy the following regularity condition Then we have the following error bounds for the discrete approximations for α = min{s − 1, 1}.

Numerical analysis of the eigenvalue problem
We now apply the Babuška-Osborn theory to derive the convergence of eigenvalues and eigenfunctions of the discrete problem (6) to those of the continuous problem (3) and estimate the order of convergence.
Using the well posedness of the continuous source problem (11), the operators T : L 2 (Ω) → H 1 0 (Ω) and S : L 2 (Ω) → L 2 0 (Ω) are well defined for any f ∈ L 2 (Ω) such that T f = u f and Sf = p f are the velocity and pressure components of the solution to problem (11).
Since the discrete source problem ( 13) is well posed, we define in the same manner the operators h and S h f = p f h are the discrete velocity and the discrete pressure approximations.Note that the operator T h is well defined in L 2 (Ω) but not in H 1 0 (Ω).Hence, we can only conclude convergence of the operators T h in L 2 (Ω) from the abstract theory.
From the a priori estimates ( 15) and ( 16) for the soure problem, we conclude which leads to the convergence of eigenvalues and eigenfunctions.From the Babuška-Osborn theory we get the following rates of convergence for the L 2 (Ω) error for the velocity component u and the L 2 (Ω) error for the pressure component p of eigenfunctions under the regularity condition ( 14) , and α = min{s − 1, 1}.From Mercier et al. [28], we conclude that the eigenvalues converge twice as fast as the eigenfunctions, i.e.
Theorem 3.2.The following estimate holds for Proof.Let λ be the eigenvalue corresponding to the eigenfunction u.Then it holds that The first two terms of the right hand side are directly estimated by ( 20) and ( 15), and the last term is estimated using (13) as follows Applying Cauchy-Schwartz inequality and ( 19), implies where α = min{s − 1, 1}.
We will now establish a relationship between the eigenvalue and the eigenfunction errors.In order to do so, we observe that the numerical scheme is consistent.
Proof.The result follows from the consistency of the discontinuous Galerkin finite element method for the source problem [32,Lemma 7.5].

then the Rayleigh quotient satisfies the following identity
Moreover, from consistency we get Next we write the following identity Subtracting ( 22) from ( 21), we obtain Dividing by (u h , u h ) on both sides in the above equation ends the proof.

A posteriori error analysis
In this section, we present a residual-based a posteriori error estimator for the Stokes eigenvalue problem.
and the edge residual estimator η E K by where I denotes the 2×2 identity matrix.Next, we introduce the estimator η J K , which measures the jump of the approximate solution u h , The local error indicator, which is the sum of the above three terms, is defined as Finally, we introduce the (global) a posteriori error estimator

Reliability
First we define the discontinuous As in [17,22], we define Hence, we decompose the DG velocity approximation uniquely into where u c h ∈ V c h and u r h ∈ V ⊥ h .Using the triangle inequality, we can write and from [17, Proposition 4.1] we get the upper bound for the second term Note that the DG form a h (u, v) is not well defined for functions u, v which belong to H 1 0 (Ω).One can overcome this difficulty by the use of a suitable lifting operator, cf.[10,22].Here, we discuss a different approach where the DG form a h (•, •) is split into several parts, Applying Cauchy-Schwarz inequality, implies Using a trace estimate together with a discrete inverse inequality leads for an edge E ∈ E h , with

Thus we have
Let Π h : H 1 0 → V c h denote the Scott-Zhang interpolation operator [30], which is stable ∇(Π h v) 0 ∇v 0 and satisfies the following interpolation property for any v ∈ H 1 0 (Ω).
Proof.Using integration by parts on each element K ∈ T h , we have Cauchy-Schwarz inequality and ( 30), lead to |||v|||.
Since (v − v c h )| ∂Ω = 0 we can rewrite T 2 in terms of a sum over interior edges Again, applying Cauchy-Schwarz inequality and (30), imply Combining the above estimates, proves the desired result.
6), then we have the following upper bound for the conforming velocity and pressure errors Proof.Using Lemma 4.1, there exists a pair (v, q) ∈ H and From (3), we obtain Applying the fact (q, ∇ • u h ) = 0, implies where Using Lemma 4.3, we have Cauchy-Schwarz inequality and (29) show Using Lemma 4.2 for the bound of T 3 , we have Cauchy-Schwarz and Poincare inequality lead to Combining the above with the estimate |||(v, q)||| |||(u − u c h , p − p h )||| yields the desired result.
Theorem 4.5.Let (u, p, λ) ∈ H 1 0 (Ω) × L 2 0 (Ω) × R + be the solution of the Stokes eigenvalue problem (3) and (u h , p h , λ h ) ∈ V h × Q h the H div -DG approximation obtained by (6).Let η h be the a posteriori error estimator in (23).Then we obtain the following a posteriori error bound where the hidden constant is independent of the viscosity ν and the sufficiently large penalty parameter γ ≥ 1.
Proof.The proof follows directly from a combination of Lemma 4.4 and (29).
The consistency term can further be estimated as where From the estimates in [29, Section 8], we can conclude that |C h (u − u h , u − u h )| is of the same order as |||(u − u h , p − p h )||| 2 .Hence C * can be bounded from above by a uniform constant.The assertion then follows from a combination of the above with Theorems 3.4 & 4.5.

Efficiency
This section is devoted to prove an efficiency bound for η.To prove the results, we use the bubble function technique which was introduced in [33,34].
Let K be an element of T h .We consider the standard element bubble function b K on K. Let v h be any vector valued polynomial function on K, then the following results hold from [1,22,33], Lemma 4.7.For u h ∈ V h , it holds that Summing over all K ∈ T h , we have Proof.Define the functions R and W locally for any K ∈ T h by From (32) we have Note that λu + ν∆u − ∇p = 0. Subtracting this from the last term, using integration by parts and W | ∂K = 0, we obtain Applying Cauchy-Schwarz inequality, implies From (32) we get Hence, dividing (33) by η R K and taking the square-root of the sum of the squares over all K ∈ T h ends the proof.
Let E be an interior edge which is shared by two elements K 1 and K 2 .Let b E denote the standard polynomial edge bubble function for E with support in ω E = {K 1 , K 2 }.In case of a regular edge E, we choose K = K 2 .When one vertex of E is a hanging node, then we choose K 1 such that E is an entire edge of K 1 and define K ⊂ K 2 as the largest rectangle contained in K 2 such that E is one of the entire edges of K. We then set ω E = {K, K}.
If σ is a vector-valued polynomial function on E, then Moreover we can define an extension σ b ∈ H 1 0 ( ω E ) such that σ b | E = b E σ and from [33,1,22] we have Proof.Let for any interior edge E ∈ E i h the functions R and Λ be such that Using (34) and [[(pI − ∇u)n]]| E = 0 we get Using Green's formula over each of the two element of ωE , gives Using λu + ν∆u − ∇p = 0, we obtain Using Cauchy-Schwarz inequality, shape-regularity of the mesh, and (35) yields as well as Combining the above estimates T 1 , T 2 and T 3 , dividing (36) by ν −1/2 h 1/2 E R 0,E and summing over all interior edges of all K ∈ T h the desired result is proven by the finite overlap of the patches ω E and Lemma 4.8.
Theorem 4.10.Let (u, p, λ) ∈ H 1 0 (Ω) × L 2 0 (Ω) × R + be the solution of the Stokes eigenvalue problem (3) and (u h , p h , λ h ) ∈ V h × Q h × R + the H div -DG approximation obtained by (6).Then the a posteriori error estimator η h is efficient in the sense that where Proof.The statement follows from a combination of Lemma 4.7-4.9.
Corollary 4.11.If u ∈ H 2 (T h ) and p ∈ H 1 (T h ), then the eigenvalue error satisfies Proof.Combining Theorems 3.4 & 4.10 implies Let C # be defined as From the estimates in [29, Section 8] and the eigenvalue estimate (20), we can conclude that ν Hence C # can be bounded from above by a uniform constant.Then it holds that This completes the proof.

Numerical experiments
This section is devoted to several numerical experiments on one convex and two non-convex domains.The experiments verify reliability and efficiency of the proposed a posteriori error estimator of Section 4 for the eigenvalue error of the smallest (simple) eigenvalue and up to polynomial degree 3. We employ the standard adaptive finite element loop with the steps solve, estimate, mark and refine.To solve the algebraic eigenvalue problem we use the ARPACK library [24] in combination with a direct solver.We mark elements of the mesh for refinement on the level in a minimal set M using the bulk marking strategy [11] with bulk parameter θ = 1/2, i.e.M is the minimal set such that θ K∈T η 2 K ≤ K∈M η 2 K .The mesh is refined with one level irregular nodes.The implementation of the method is done in the software library amandus [20], which is based on the dealii finite element library [5].
In all experiments we consider the viscosity ν = 1 and chose the penalty parameter γ = k(k + 1)/2 for k-th order RT k × Q k finite element pairs, k = 1, 2, 3. Since the eigenvalues of the Stokes problem are related to the eigenvalues of the buckling eigenvalue problem of clamped plates via the stream function formulation, we can use reference values for the eigenvalues from [6,7,31].

Square domain
In this example, we consider the square domain Ω = (0, 1) 2 .The reference value for the first eigenvalue reads λ = 52.344691168[6,7,31].The streamline plot of the discrete eigenfunction u and the plot of the discrete pressure p on a uniform mesh for k = 1 are displayed in Figures 1(a) and 1(b), respectively.In Figure 2, we observe that both uniform and adaptive mesh refinement leads to optimal orders of convergence O(N −k ) for the eigenvalue error |λ − λ |.This is due to the fact that the domain is convex and the first eigenfunction is smooth enough.Note that for uniform meshes O(N −k ) ≈ O(h 2k ), for N = dim(V h × Q h ).We observe that the convergence graphs for uniform and adaptive mesh refinement overlap each other for both the eigenvalue errors |λ − λ | as well the a posteriori error estimators η 2 .Moreover, we confirm that the a posteriori error estimator η 2 is numerically reliable and efficient.

L-shaped domain
In the second example, we take the non-convex L-shaped domain Ω = (−1, 1) 2 \ (0, 1) 2 with a re-entrant corner at the origin, which allows for singular eigenfunctions.To compute the error of the first eigenvalue, we take λ = 32.13269465 as a reference value.Figures 3(a) and 3(b) show the computed velocity and discrete pressure as a streamline plot on a uniform mesh computed with k = 1.The exponent for the singular function at the re-entrant corner is known to be α ≈ 0.544483736782464.Hence, in Figure 4 we observe suboptimal convergence of O(N −0.544 ) for the eigenvalue error even for k = 2. Adaptive mesh refinement however, achieves optimal convergence O(N −k ) of the eigenvalue error for k = 1, 2, 3.The a posteriori error estimator η 2 shows to be reliable and efficient in all experiments.Observe that the eigenvalue error obtained with k = 3 on adaptively refined meshes is about 6 orders of magnitude smaller than that for uniform mesh refinement.This demonstrates the importance of mesh adaptivity, in particular for high order methods.Figures 5(a) and 5(b) show some adaptively refined meshes for k = 1, 2, 3, which show strong refinement towards the origin.

Slit domain
In the last example, let Ω = (−1, 1) 2 \ ({0} × (−1, 0)) be the slit domain.To compute the error of the first eigenvalue, we take λ = 29.9168629as a reference value.The discrete velocity eigenfunction u and discrete pressure p are displayed in Figures 6(a) and 6(b) as a streamline plot on a uniform mesh for k = 1.In Figure 7, we observe suboptimal convergence of O(N −1/2 ) for the eigenvalue error on uniform meshes, but optimal convergence of O(N −k ) for k = 1, 2, 3, for adaptively refined meshes.Moreover, the a posteriori error estimator proves to be numerically reliable and efficient.Note that we had to stop the third order method on adaptively refined 2 and τ = (τ ij ) 2×2 , then

Figure 1 :
Figure 1: Streamline plot of the discrete eigenfunction u (a), and plot of the discrete pressure p (b).

Figure 2 :
Figure 2: Convergence history of |λ − λ | and η 2 on uniformly and adaptively refined meshes for the square domain.

Figure 3 :
Figure 3: Streamline plot of the discrete eigenfunction u (a) , and plot of the of discrete pressure p (b).

Figure 4 :
Figure 4: Convergence history of |λ − λ | and η 2 on uniformly and adaptively refined meshes for the L-shaped domain.

Figure 6 :
Figure 6: Streamline plot of the discrete eigenfunction u (a) , and plot of the discrete pressure p (b).

Figure 7 :
Figure 7: Convergence history of |λ − λ | and η 2 on uniformly and adaptively refined meshes for the slit domain.
17, Lemma 4.3], [22, Section 2.3].We include the proof for the H div -DG formulation of the Stokes problem for completeness.