A stabilizer free WG method for the Stokes equations with order two superconvergence on polytopal mesh

1.   Introduction

A stabilizing/penalty term is often used in finite element methods with discontinuous approximations to enforce connection of discontinuous functions across element boundaries. Development of stabilizer free discontinuous finite element method is desirable since it simplifies finite element formulation and reduces programming complexity. The stabilizer free WG method and the stabilizer DG method on polytopal mesh were first introduced in [11,12] for second order elliptic problems. The main idea in [11,12] is to raise the degree of polynomials used to compute weak gradient $\nabla_w$. In [11,12], gradient is approximated by a polynomial of order $j = k+n-1$ where $n$ is the number of sides of polygonal element. This result has been improved in [1,2] by reducing the degree of polynomial $j$. Recently, new stabilizer free WG methods have been developed in [13,14] for second order elliptic equations on polytopal mesh, which have superconvergence. Wachspress coordinates are used to approximate $\nabla_w$ in [6,7] for solving the Stokes equations on polytopal mesh. Wachspress coordinates are usually rational functions, instead of polynomials. The WG methods in [6,7] are limited to the lowest order WG elements.

In general, discontinuous finite element methods tend to have complex formulations which are often necessary to enforce weak continuity of discontinuous solutions across element boundaries. Most of discontinuous finite element methods have one or more stabilizing terms to guarantee stability and convergence of the methods. The stabilizer free WG method has a super clean finite element formulations as (4)-(5) compared to the weak Galerkin method and other DG methods, which reduces programming complexity.

One obvious disadvantage of discontinuous finite element methods is their rather complex formulations which are often necessary to enforce weak continuity of discontinuous solutions across element boundaries. Most of discontinuous finite element methods have one or more stabilizing terms to guarantee stability and convergence of the methods. Existing of stabilizing terms further complicates formulations.

In this paper, we introduce a new stabilizer free WG method of any order to solve the Stokes problem: find unknown functions ${\bf u}$ and $p$ such that

where $\Omega$ is a polygonal or polyhedral domain in $\mathbb{R}^d\; (d = 2,3)$. Our new WG method has the following formulations without any stabilizers: seek $( {\bf u}_h,p_h)\in V_h\times W_h$ satisfying the following for all $( {\bf v},w)\in V_h\times W_h$,

Here $\nabla_w$ and $\nabla_w\cdot$ are weak gradient and weak divergence, respectively. In addition, we have proved that the WG approximations have the convergence rates two order higher than the optimal-order for velocity in both an energy norm and the $L^2$ norm and the optimal convergence rate for pressure in the $L^2$ norm. Extensive numerical examples are tested for the new WG elements of different degrees $k$ in both two and three dimensional spaces.

Comparing to other weak Galerkin finite element methods, the WG methods without stabilizers leads to a 2-order superconvergent solution with a carefully chose weak Galerkin. So far all weak Galerkin finite element methods with stabilizers do not have superconvergence in both an energy norm and the $L^2$ norm. A order one superconvergence in an energy norm only is derived in [4,3]. The new method reduces the computation cost greatly of the other WG methods. The price to pay for removing the stabilizers from the WG methods is that the weak gradient, an intermediate variable, need to be approximated by one degree higher and piecewise polynomials.

2.   Preliminary

Let ${\mathcal T}_h$ be a partition of the domain $\Omega$ consisting of polygons in two dimension or polyhedra in three dimension satisfying a set of conditions specified in [10]. Denote by ${\mathcal E}_h$ the set of all edges or flat faces in ${\mathcal T}_h$, and let ${\mathcal E}_h^0 = {\mathcal E}_h\backslash\partial\Omega$ be the set of all interior edges or flat faces. For every element $T\in {{\mathcal T}}_h$, we denote by $h_T$ its diameter and mesh size $h = \max_{T\in{{\mathcal T}}_h} h_T$ for ${\mathcal T}_h$. Let $P_k(T)$ consist all the polynomials on $T$ with degree no greater than $k$.

For $k\ge 0$ and given ${{\mathcal T}}_h$, define two finite element spaces for velocity

and for pressure

Let $V_h^0$ be a subspace of $V_h$ consisting of functions with vanishing boundary value.

The space $H(\operatorname{div};\Omega)$ is defined as

For any $T\in{{\mathcal T}}_h$, it can be divided in to a set of disjoint triangles $T_i$ with $T = \cup T_i$. Then we define a space $\Lambda_h(T)$ for the approximation of weak gradient on each element $T$ as

For a function ${\bf v}\in V_h$, its weak gradient $\nabla_w {\bf v}$ is a piecewise polynomial satisfying $\nabla_w {\bf v}|_T \in \Lambda_k(T)$ and the following equation,

For a function ${\bf v}\in V_h$, its weak divergence $\nabla_w\cdot {\bf v}$ is a piecewise polynomial satisfying $\nabla_w\cdot {\bf v}|_T \in P_{k+1}(T)$ and the following equation,

The proof of the following lemma can be found in [14].

Lemma 2.1. For $\tau\in [H(\operatorname{div};\Omega)]^d$, there exists a projection $\Pi_h$ with $\Pi_h\tau\in [H(\operatorname{div};\Omega)]^d$ satisfying $\Pi_h\tau|_T\in \Lambda_k(T)$ and the followings

3.   Finite element method and its well posedness

We start this section by introducing the following WG finite element scheme without stabilizers.

Weak Galerkin Algorithm 3.1. A numerical approximation for (1)-(3) is finding $( {\bf u}_h,p_h)\in V_h^0\times W_h$ such that for all $( {\bf v},w)\in V_h^0\times W_h$,

Let $Q_0$ and $Q_b$ be the two element-wise defined $L^2$ projections onto $[P_k(T)]^d$ and $[P_{k+1}(e)]^d$ with $e\subset\partial T$ on $T$ respectively. Define $Q_h {\bf u} = \{Q_0 {\bf u},Q_b {\bf u}\}\in V_h$ for the true solution ${\bf u}$. Let ${{\mathbb Q}}_h$ be the element-wise defined $L^2$ projection onto $\Lambda_k(T)$ on each element $T$. Finally denote by $\mathcal{Q}_h$ the element-wise defined $L^2$ projection onto $P_{k+1}(T)$ on each element $T$.

Lemma 3.1. Let $\boldsymbol\phi\in [H_0^1(\Omega)]^d$, then on $T\in{{\mathcal T}}_h$

Proof. Using (8) and integration by parts, we have that for any $\tau\in \Lambda_k(T)$

which implies the identity (15).

Using (9) and integration by parts, we have that for any $w\in P_{k+1}(T)$

which proves (15).

For any function $\varphi\in H^1(T)$, the following trace inequality holds true (see [10] for details):

We introduce two semi-norms ${{|||}} {\bf v}{{|||}}$ and $\| {\bf v}\|_{1,h}$ for any ${\bf v}\in V_h$ as follows:

It is easy to see that $\| {\bf v}\|_{1,h}$ defines a norm in $V_h^0$. Next we will show that ${{|||}} \cdot {{|||}}$ also defines a norm in $V_h^0$ by proving the equivalence of ${{|||}}\cdot{{|||}}$ and $\|\cdot\|_{1,h}$ in $V_h$.

The following norm equivalence has been proved in [14] for each component of ${\bf v}$,

Unlike the traditional finite elements [5,8,9,16,17,18,19,20,21], the inf-sup condition for the weak Galerkin finite element is easily satisfied due to the large velocity space with independent element boundary degrees of freedom.

Lemma 3.2. There exists a positive constant $\beta$ independent of $h$ such that for all $\rho\in W_h$,

Proof. For any given $\rho\in W_h\subset L_0^2(\Omega)$, there exists a function $\tilde {\bf v}\in [H_0^1(\Omega)]^d$ such that

where $C>0$ is a constant independent of $h$. Let ${\bf v} = Q_h\tilde{ {\bf v}} = \{Q_0\tilde{ {\bf v}},Q_b\tilde{ {\bf v}}\}\in V_h$. It follows from (20), (17) and $\tilde{ {\bf v}}\in [H_0^1(\Omega)]^d$,

which implies

It follows from (9) that

Using (24), (23) and (22), we have

for a positive constant $\beta$. This completes the proof of the lemma.

Lemma 3.3. The weak Galerkin method (13)-(14) has a unique solution.

Proof. It suffices to show that zero is the only solution of (13)-(14) if ${{\bf f}} = 0$. To this end, let ${{\bf f}} = 0$ and take ${\bf v} = {\bf u}_h$ in (13) and $w = p_h$ in (14). By adding the two resulting equations, we obtain

which implies that $\nabla_w {\bf u}_h = 0$ on each element $T$. By (20), we have $\| {\bf u}_h\|_{1,h} = 0$ which implies that ${\bf u}_h = 0$.

Since ${\bf u}_h = 0$ and ${{\bf f}} = 0$, the equation (13) becomes $(\nabla\cdot {\bf v},\ p_h) = 0$ for any ${\bf v}\in V_h$. Then the inf-sup condition (21) implies $p_h = 0$. We have proved the lemma.

4.   Error equations

In this section, we derive the equations that the errors satisfy. Let ${\bf e}_h = Q_h {\bf u}- {\bf u}_h$ and $\varepsilon_h = \mathcal{Q}_hp-p_h$.

Lemma 4.1. The following error equations hold true for any $( {\bf v},w)\in V_h^0\times W_h$,

where

Proof. First, we test (1) by ${\bf v}_0$ with ${\bf v} = \{ {\bf v}_0, {\bf v}_b\}\in V_h^0$ to obtain

It follows from (11) and (15)

Using integration by parts and the fact ${{\langle}} p, {\bf v}_b\cdot{{\bf n}}{{\rangle}}_{\partial{{\mathcal T}}_h} = 0$, we have

which implies

Substituting (30) and (31) into (29) gives

The difference of (32) and (13) implies

Testing equation (2) by $w\in W_h$ and using (16) give

The difference of (34) and (14) implies (26). We have proved the lemma.

5.   Error estimates in energy norm

In this section, we establish order two superconvergence for the velocity approximation ${\bf u}_h$ in ${{|||}}\cdot{{|||}}$ norm and optimal order error estimate for the pressure approximation $p_h$ in the standard $L^2$ norm.

Lemma 5.1. Let ${\bf u}\in [H^{k+3}(\Omega)]^d$ and $p\in H^{k+2}(\Omega)$ and ${\bf v}\in V_h$. Then, the following estimates hold true

Proof. Using the Cauchy-Schwarz inequality and the definitions of ${{\mathbb Q}}_h$ and $\Pi_h$, we have

It follows from (17) and (20)

We have proved the lemma.

Theorem 5.2. Let $( {\bf u}_h,p_h)\in V_h^0\times W_h$ be the solution of (13)-(14). Then, we have

Proof. By letting ${\bf v} = {\bf e}_h$ in (25) and $w = \varepsilon_h$ in (26) and using the equation (26), we have

It then follows from (35) and (36) that

which implies (37). To estimate $\|\varepsilon_h\|$, we have from (25) that

Using (40), (35) and (36), we arrive at

Combining the above estimate with the inf-sup condition (21) gives

which yields the desired estimate (38).

6.   Error estimates in L2 norm

In this section, order two superconvergence for velocity in the $L^2$ norm is obtained by duality argument. Recall that ${\bf e}_h = \{ {\bf e}_0, {\bf e}_b\} = Q_h {\bf u}- {\bf u}_h$ and $\epsilon_h = \mathcal{Q}_hp-p_h$. Consider the dual problem: seeking $( {\boldsymbol\psi},\xi)$ satisfying

Assume that the dual problem (41)-(43) satisfy the following regularity assumption:

We need the following lemma first.

Lemma 6.1. For any ${\bf v}\in V_h^0$ and $w\in W_h$, the following equations hold true,

where

Proof. Testing (41) by ${\bf v}_0$ with ${\bf v} = \{ {\bf v}_0, {\bf v}_b\}\in V_h^0$ gives

It follows from integration by parts and the fact $\langle\nabla {\boldsymbol\psi}\cdot{{\bf n}}, {\bf v}_b\rangle_{\partial{{\mathcal T}}_h} = 0$

By integration by parts, (8) and (15)

Combining (48) and (49) gives

Similar to the derivation of (31), we obtain

Combining (50) and (51) with (47) yields (45). Testing equation (42) by $w\in W_h$ and using (16) give

which implies (46) and we have proved the lemma.

By the same argument as (50), (25) has another form as

Theorem 6.2. Let $( {\bf u}_h,p_h)\in V_h^0\times W_h$ be the solution of (13)-(14). Assume that (44) holds true. Then, we have

Proof. Letting ${\bf v} = {\bf e}_h$ in (45) yields

Using the fact $(Q_h\xi,\;\nabla_w\cdot {\bf e}_h) = 0$, (55) becomes

With ${\bf v} = Q_h {\boldsymbol\psi}$, (53) becomes

Using (46), we have $(\epsilon_h,\;\nabla_w\cdot Q_h {\boldsymbol\psi}) = 0$. Then (57) becomes

Combining (56) and (58), we have

Using the Cauchy-Schwarz inequality, the trace inequality (17) and the definition of ${{\mathbb Q}}_h$, we arrive at

Similarly, we have

It follows from the Cauchy-Schwarz inequality, the trace inequality (20) and (37),

Similarly,

Combining all the estimates above with (59) yields

The estimate (54) follows from the above inequality and the regularity assumption (44). We have completed the proof.

7.   Numerical experiments

7.1.   Example 1

We solve the following 2D stationary Stokes equations with domain $\Omega = (0,1)^2$:

The exact solution is

In this example, we use quadrilateral grids shown in Figure 1.

We test the newly constructed, two-order superconvergent weak Galerkin finite elements. We compare the results with that of the standard optimal order convergent (no superconvergence) weak Galerkin finite elements. The two-order superconvergent weak Galerkin finite element spaces are

The standard weak Galerkin finite element spaces are

In Table 1, we list the errors and the orders of convergence, for the two types $P_1$ weak Galerkin finite elements. We can see that two-order superconvergence is achieved for the velocity in $L^2$-norm and $H^1$-like norm, for the new weak Galerkin finite element (65). The pressure converges at the optimal order in this case. But the standard weak Galerkin finite element converges only at the optimal order, with all errors 1000 times bigger than that of the two-order superconvergence weak Galerkin finite element (65).

To see the superconvergence phenomenon, in Figure 2, we plot the first component of ${\bf u}_h$ of the $P_1^2$-$P_2^2$-$P_2$ WG finite element (65) solution for (64) on the fifth grid of Figure 1 on the top. In middle of Figure 2, we plot its error. At the bottom of Figure 2, we plot the error of the $P_1^2$-$P_1^2$-$P_1$ WG finite element (66) solution $( {\bf u}_1)_h$ on the fifth grid. We note that both solutions are $P_1$ polynomials, but the error of latter is 1000 times bigger.

In Figure 3, we plot the second component of ${\bf u}_h$ of the $P_1^2$-$P_2^2$-$P_2$ WG finite element (65) solution for (64) on the fifth grid of Figure 1 on the top. In middle of Figure 3, we plot its error. At the bottom of Figure 3, we plot the error of the $P_1^2$-$P_1^2$-$P_1$ WG finite element (66) solution $( {\bf u}_1)_h$ on the fifth grid. Correspondingly, in Figure 4, we plot the numerical solution for pressure, and the two errors.

In Table 2, we compute the two-order superconvergent $P_2$ weak Galerkin finite element solutions and the standard $P_2$ weak Galerkin finite element solutions. We can see the former has two orders higher convergent rate than the latter. On same grids and with same $P_2$ polynomials for ${\bf u}_0$, the error of former is 1000 times smaller than the latter. In Table 3, both types of $P_3$ weak Galerkin finite element solutions are listed. They verify the theory.

7.2.   Example 2

We solve the problem (64) again, on polygonal grids, consisting of quadrilaterals, pentagons and hexagons, shown in Figure 5. We intentionally perturb the grids to show that the two-order superconvergence is grid independent, i.e., not necessarily on uniform grids. In Table 4, we list the errors and the orders of convergence. The computational results match the theoretic order of convergence, in all cases.

7.3.   Example 3

We compute the following 2D driven cavity flow on domain $\Omega = (0,1)^2$:

The solution is not in $H^1_0(\Omega)$. But the discontinuous ${\bf u}_b$ function can handle the jump boundary condition well. Additionally, the discontinuous ${\bf u}_b$ function can be defined well on grids with hanging nodes, i.e., non-compatible grids (see Figure 6). In this computation, we use the graded grid shown in Figure 6.

The computed velocity field is plotted in Figure 7.

The extra fine grid near the two top corners gives an accurate solution at these singularity points. The zoom-in plots at the two corners are displayed in Figures 8-9.

Because we use smaller square elements at the two bottom corners, we can compute the secondary flow at these two corners well. They are plotted in Figures 10-11.

We zoom-in further the solution at the lower-left corner. We can see in Figure 12 there is another secondary flow generated by the secondary flow in Figure 10. Nevertheless we can see such a computed flow is not mass conservative. We need to use $H(\hbox{{div}})$ finite elements [15] or divergence-free $H^1$ finite elements [5,8,9,16,17,18,19,20,21] to get a mass conservative solution.

7.4.   Example 4

We compute a 3D problem (1)-(3) with $\Omega = (0,1)^3$. The source term $\bf$ and the boundary value ${{\bf g}}$ are chosen so that the exact solution is

We use tetrahedral meshes shown in Figure 13. The results of the 3D $P_k$-$P_{k+1}$ weak Galerkin finite element methods are listed in Table 5. The results show that the method is stable and is of two-order superconvergence (for velocity).