Estimation and numerical validation of inf-sup constant for bilinear form $(p,\mathrm{div}\ \bm{u})$

Numerical schemes for mixed finite element methods should result in stable bounded solutions. The stability of finite element discretization is governed by the inf-sup criterion for bilinear forms. The objective of this paper is to derive and validate the inf-sup stability constant for the bilinear form of the divergence term. This is of importance because it often appears as pressure constraint term in discretization of fluid mechanics equations. For given $\bm{u}\in H\left(\mathrm{div};\Omega\right)$, $p\in L^{2}\left(\Omega\right)$, the bilinear form is give by

The inf-sup constant $\beta$ for (1) is given by

Quite some work has been done to approximate the value of inf-sup constant see for eg. [2, 3, 4]. In this work we measure the norm of the velocity field in the space orthogonal to kernel of the divergence operator and prove that the value of the inf-sup constant $\beta=1.0$. For numerical tests we use finite dimensional spaces defined in [1]. The validation study is performed on three different domains: i) unit square $\Omega=[0,1]^{2}$; ii) square $\Omega=[0,2]^{2}$, for varying mesh sizes, $h$, and polynomial degrees $N=1,2,3$. The numerical values are in agreement with the derived values of $\beta$, with maximum errors shown of the order of $10^{-6}$.

Let $\Omega\subset\mathbb{R}^{2}$ be an open, bounded domain. We will use the finite dimensional spaces and the divergence operator defined in [1, §3]: $D\left(\Omega\right)\subset H\left(\mathrm{div};\Omega\right)$, $\widetilde{S}\left(\Omega\right)\subset L^{2}\left(\Omega\right)$, $S\left(\Omega\right)\subset L^{2}\left(\Omega\right)$, and $\mathbb{E}^{2,1}$ the discrete representation of the divergence operator. Let $K=\mathrm{Ker}\ \mathbb{E}^{2,1}$, $H=\mathrm{Ker}\left({\mathbb{E}^{2,1}}^{\intercal}\right)$. The discrete inf-sup condition is then given by

where $\mathcal{N}^{x}$ denotes the degrees of freedom in our finite element, i.e. the vector of expansion coefficients. Here $p\in\widetilde{S}(\Omega)\subset L^{2}(\Omega)$ and $\bm{u}\in D(\Omega)\subset H(\mathrm{div};\Omega)$ and the vectors $\widetilde{\mathcal{N}}^{0}\left(p\right)$ and $\mathcal{N}^{1}\left(\bm{u}\right)$ are the expansion coefficients. The norm of $p\in\widetilde{S}(\Omega)$ is

the norm of $\bm{u}\in H\left(\mathrm{div};\Omega\right)$ is

and the norm of $\bm{u}\in K^{\perp}\subset D(\Omega)$ is

In the continuous case, using Cauchy Schwartz inequality, we have, for all $p\in L^{2}\left(\Omega\right)$ and $\bm{u}\in K^{\perp}$

If this inequality holds for all $p$ and $\bm{u}$ it should also hold when we take the infimum over $\widetilde{S}(\Omega)$ and the supremum over $D(\Omega)$, from which we conclude that $\beta\leq 1$.

Now, for an arbitrary vector field $\bm{u}^{*}$, we have

If we now take for $\bm{u}^{*}$, the vector field with expansion coefficients which satisfy $\mathbb{M}^{(2)}\mathbb{E}^{2,1}\mathcal{N}^{1}\left(\bm{u}^{*}\right)=\widetilde{\mathcal{N}}^{0}\left(p\right)$, then $\|p\|_{L^{2}\left(\Omega\right)}=\|\bm{u}^{*}\|_{H(\mathrm{div};\Omega)}$ and the numerator $\mathcal{N}^{1}\left(\bm{u}\right)^{\intercal}\mathbb{E}^{\intercal}\widetilde{\mathcal{N}}^{0}\left(p\right)=\|p\|_{L^{2}\left(\Omega\right)}^{2}=\|\bm{u}^{*}\|_{H(\mathrm{div};\Omega)}\|p\|_{L^{2}\left(\Omega\right)}$. If we insert these estimates in (7) we have

which shows that $\beta\geq 1$. From (6) and (8) we conclude that $\beta=1$ in the discrete setting. This value is independent of the mesh size or polynomial degree, so this value is also the inf-sup constant for $h\rightarrow 0$. This value is also independent of the domain $\Omega$.

In this section we will follow [5, §3.4.3] to evaluate the inf-sup constant. Let $S_{x}$ and $S_{y}$ be the two symmetric and (semi-)positive definite matrices, such that

We can write the norms in (4) and (5) as

where $||\cdot||_{E}$ is the Eucledian vector norm.

Let $M=S_{y}^{-1}\ \mathbb{E}^{2,1}\ S_{x}^{-1}$ and its singular value decomposition be given by

Now we can write $\mathbb{E}^{2,1}$ as

Also, let

Now, if we substitute $\bm{u}$, $p$ from (15) and $\mathbb{E}^{2,1}$ from (14) in the RHS term of (3), we get

where in the last step we used [5, Prop 3.4.3] which states that there exists a positive constant ${\beta}_{h}$ that is equivalent to the smallest positive singular value of the matrix, $M=S_{y}^{-1}E^{2,1}S_{x}^{-1}$.

In Figure 1, and Figure 2, we plot the value of the inf-sup term obtained using (16) for domain $\Omega=[0,1]^{2}$, and $\Omega=[0,2]^{2}$, respectively. On the y-axis we have $\beta_{h}$, and on the x-axis we have the length of the element, $h=1/K$ for $K=1\;,2\;,4\;,8\;,16\;,32\;,64$. The plots on the right side in the figure are a zoom in view of the plots on the left side. The numeric values for plots in Figure 1 and Figure 2 are given in Table 1 and Table 2 respectively. In both the figures and the plots we observe that the numeric values of inf-sup term are all very close to $1.0$ and exact upto at least six decimal places, which is in agreement with theoretical derivation.

In this paper, we derive a theoretical estimate for the discrete inf-sup formulation and validate the value of the constant using finite dimensional spaces defined in [1]. The theoretical proof of the inf-sup term becomes straight forward when we use the appropriate norm on $K^{\perp}$ space, see (5). We evaluate the constant for two different test cases, i) a unit square domain $\Omega=[0,1]^{2}$, ii) a square domain $\Omega=[0,2]^{2}$. It is shown that for all the cases the numerical value is in agreement with the theoretical value.