Well-posedness and H(div)-conforming finite element approximation of a linearised model for inviscid incompressible flow

TL;DR

This paper establishes the well-posedness of a linearised inviscid flow model and develops H(div)-conforming finite element methods, providing error estimates for velocity and pressure approximations.

Contribution

It introduces a regularisation approach for the model and proves error bounds for finite element discretisations, advancing numerical analysis of inviscid flow models.

Findings

Existence and uniqueness of weak solutions for smooth domains.

Error estimate of order O(h^{k+1/2}) for velocity in L2-norm.

Error estimates for pressure in L2-norm.

Abstract

We consider a linearised model of incompressible inviscid flow. Using a regularisation based on the Hodge Laplacian we prove existence and uniqueness of weak solutions for smooth domains. The model problem is then discretised using H(div)-conforming finite element methods, for which we prove error estimates for the velocity approximation in the $L^2$-norm of order $O(h^{k+\frac12})$. We also prove error estimates for the pressure error in the $L^2$-norm.

Tables

Table 1. Table 1. Errors for the BDM1/P0 element. σ = 100 𝜎 100 \sigma=100 n = 1 𝑛 1 n=1 .

$h$	${\Vert 𝒖-{𝒖}_h\Vert }_{L^2\left(\Omega \right)}$	${\Vert p-p_h\Vert }_{L^2\left(\Omega \right)}$	CPU
$1/10$	$0.011(-)$	$0.15(-)$	$0.073s$
$1/20$	$0.0030\left(1.9\right)$	$0.074\left(1.0\right)$	$0.47s$
$1/40$	$0.00087\left(1.8\right)$	$0.037\left(1.0\right)$	$4.7s$
$1/80$	$0.00031\left(1.5\right)$	$0.019\left(1.0\right)$	$62.9s$

Table 2. Table 2. Errors for the RT1/P1dc element. σ = 100 𝜎 100 \sigma=100 n = 1 𝑛 1 n=1 .

$h$	${\Vert 𝒖-{𝒖}_h\Vert }_{L^2\left(\Omega \right)}$	${\Vert p-p_h\Vert }_{L^2\left(\Omega \right)}$	CPU
$1/10$	$0.011(-)$	$0.026(-)$	$0.17s$
$1/20$	$0.0030\left(1.9\right)$	$0.0060\left(2.1\right)$	$1.2s$
$1/40$	$0.00087\left(1.8\right)$	$0.0018\left(1.7\right)$	$12s$
$1/80$	$0.00031\left(1.5\right)$	$0.00073\left(1.3\right)$	$165s$

Table 3. Table 3. Errors for the BDM1/P0 element (columns 2 and 3) and RT1/P1dc element (columns 4 and 5), h = 1 / 40 ℎ 1 40 h=1/40 , σ = 100 𝜎 100 \sigma=100 , varying n 𝑛 n .

$n$	BDM ${\Vert 𝒖-{𝒖}_h\Vert }_{L^2\left(\Omega \right)}$	BDM ${\Vert p-p_h\Vert }_{L^2\left(\Omega \right)}$	RT ${\Vert 𝒖-{𝒖}_h\Vert }_{L^2\left(\Omega \right)}$	RT ${\Vert p-p_h\Vert }_{L^2\left(\Omega \right)}$
$1$	$0.00087$	$0.037$	$0.00087$	$0.0018$
$2$	$0.0048$	$0.074$	$0.0048$	$0.0058$
$4$	$0.031$	$0.14$	$0.031$	$0.026$
$8$	$0.21$	$0.34$	$0.21$	$0.18$

Table 4. Table 4. Errors for the BDM1/P0 element (columns 2 and 3) and RT1/P1dc element (columns 4 and 5), h = 1 / 40 ℎ 1 40 h=1/40 , n = 1 𝑛 1 n=1 , varying σ 𝜎 \sigma .

$\sigma$	BDM ${\Vert 𝒖-{𝒖}_h\Vert }_{L^2\left(\Omega \right)}$	BDM ${\Vert p-p_h\Vert }_{L^2\left(\Omega \right)}$	RT ${\Vert 𝒖-{𝒖}_h\Vert }_{L^2\left(\Omega \right)}$	RT ${\Vert p-p_h\Vert }_{L^2\left(\Omega \right)}$
${10}^6$	$0.00061$	$0.037$	$0.00061$	$0.015$
$100$	$0.00087$	$0.037$	$0.00087$	$0.0018$
$50$	$0.0012$	$0.037$	$0.0012$	$0.0019$
$25$	$0.0021$	$0.037$	$0.0021$	$0.0022$
$10$	$0.0051$	$0.037$	$0.0051$	$0.0045$
$1$	$0.048$	$0.058$	$0.048$	$0.045$

Equations

$$\Vert u-u_h{\Vert }_{L^2(\Omega )}\le C{h}^{k+\frac{1}{2}}|u{|}_{H^{k+1}(\Omega )},$$

$$\mathrm{div}(\mathbf{u}\otimes \mathbf{\beta })+\sigma \mathbf{u}+\nabla p=$$

$$\mathrm{div}\mathbf{u}=$$

$$\mathbf{u}\cdot \mathbf{n}=$$

$$\tilde{\mathbf{f}}(\mathbf{x},t):=\mathbf{f}(\mathbf{x})\exp (\sigma t),\sigma \in \mathbb{R}\setminus 0.$$

$${\partial }_t\tilde{\mathbf{u}}+\mathrm{div}(\tilde{\mathbf{u}}\otimes \mathbf{\beta })+\mathrm{div}(\mathbf{\beta }\otimes \tilde{\mathbf{u}})+\nabla \tilde{p}=\tilde{\mathbf{f}}(\mathbf{x},t).$$

$$\tilde{\mathbf{u}}(\mathbf{x},t)=\mathbf{u}(\mathbf{x})\exp (\sigma t).$$

$$\sigma \mathbf{u}+\mathrm{div}(\mathbf{u}\otimes \mathbf{\beta })+\mathrm{div}(\mathbf{\beta }\otimes \mathbf{u})+\nabla p=\mathbf{f}(\mathbf{x}).$$

$${\bm{\beta}}:=\left(\begin{array}[]{c}\varphi(y)\\ 0\end{array}\right).$$

$${\bm{\beta}}:=\left(\begin{array}[]{c}\partial_{y}\varphi(x,y)\\ -\partial_{x}\varphi(x,y)\end{array}\right).$$

$$p=n^2{\pi }^2({\cos }^2(n\pi x)-{\sin }^2(n\pi y))/2.$$

$$\mathrm{div}\mathbf{v}=q\text{ in }\Omega ,$$

$$\Vert \nabla \mathbf{v}{\Vert }_{\Omega }\le C\Vert q{\Vert }_{\Omega }.$$

$$\Vert \mathbf{v}{\Vert }_{\Omega }\le C(\Vert \mathrm{div}\mathbf{v}{\Vert }_{\Omega }+\Vert \mathrm{curl}\mathbf{v}{\Vert }_{\Omega })\forall \mathbf{v}\in {\mathbf{H}}_0(\mathrm{div},\Omega )\cap \mathbf{H}(\mathrm{curl},\Omega ).$$

$$\Vert \nabla \mathbf{v}{\Vert }_{\Omega }\le K(\Vert \mathrm{div}\mathbf{v}{\Vert }_{\Omega }+\Vert \mathrm{curl}\mathbf{v}{\Vert }_{\Omega })\forall \mathbf{v}\in {\mathbf{H}}_0(\mathrm{div},\Omega )\cap \mathbf{H}(\mathrm{curl},\Omega ).$$

$$\Vert \nabla \mathbf{v}{\Vert }_{\Omega }^2\le \Vert \mathrm{div}\mathbf{v}{\Vert }_{\Omega }^2+\Vert \mathrm{curl}\mathbf{v}{\Vert }_{\Omega }^2\forall \mathbf{v}\in {\mathbf{H}}_0(\mathrm{div},\Omega )\cap \mathbf{H}(\mathrm{curl},\Omega ).$$

$$\mathrm{curl}(\mathbf{\beta }\cdot \nabla \mathbf{v})=\mathbf{\beta }\cdot \nabla (\mathrm{curl}\mathbf{v})+((\nabla \mathbf{\beta }{)}^t\times \nabla \mathbf{v}).$$

$$-(\mathbf{u},\mathbf{\beta }\cdot \nabla \mathbf{v})+(\sigma \mathbf{u},\mathbf{v})-(p,\mathrm{div}\mathbf{v})=$$

$$(\mathrm{div}\mathbf{u},q)=$$

$$-\epsilon \Delta {\mathbf{u}}_{\epsilon }+\mathrm{div}({\mathbf{u}}_{\epsilon }\otimes \mathbf{\beta })+\sigma {\mathbf{u}}_{\epsilon }+\nabla p_{\epsilon }=$$

$$\mathrm{div}{\mathbf{u}}_{\epsilon }=$$

$${\mathbf{u}}_{\epsilon }=$$

$$\epsilon (\nabla {\mathbf{u}}_{\epsilon },\nabla \mathbf{v})-({\mathbf{u}}_{\epsilon },\mathbf{\beta }\cdot \nabla \mathbf{v})+(\sigma {\mathbf{u}}_{\epsilon },\mathbf{v})-(p_{\epsilon },\mathrm{div}\mathbf{v})=$$

$$(\mathrm{div}{\mathbf{u}}_{\epsilon },q)=$$

$$\Vert p_{\epsilon }{\Vert }_{\Omega }+\Vert \sqrt{\sigma }{\mathbf{u}}_{\epsilon }{\Vert }_{\Omega }+\sqrt{\epsilon }\Vert \nabla {\mathbf{u}}_{\epsilon }{\Vert }_{\Omega }\le C\Vert \mathbf{f}{\Vert }_{\Omega },$$

$$\epsilon \Vert \nabla {\mathbf{u}}_{\epsilon }{\Vert }_{\Omega }^2+\Vert \sqrt{\sigma }{\mathbf{u}}_{\epsilon }{\Vert }_{\Omega }^2=(\mathbf{f},{\mathbf{u}}_{\epsilon }).$$

$$\epsilon \Vert \nabla {\mathbf{u}}_{\epsilon }{\Vert }_{\Omega }^2+\frac{1}{2}\Vert \sqrt{\sigma }{\mathbf{u}}_{\epsilon }{\Vert }_{\Omega }^2\le \frac{1}{2{\sigma }_0}\Vert \mathbf{f}{\Vert }_{\Omega }^2.$$

$$\|p_{\varepsilon}\|_{\Omega}\leq C\,\Big{(}\varepsilon\|\nabla\bm{u}_{\varepsilon}\|_{\Omega}+\|{\bm{\beta}}\|_{\infty,\Omega}\|\bm{u}_{\varepsilon}\|_{\Omega}+\|\sqrt{\sigma}\|_{\infty,\Omega}\,\|\sqrt{\sigma}\,\bm{u}_{\varepsilon}\|_{\Omega}+\|\bm{f}\|_{\Omega}\Big{)}\,,$$

$$\epsilon \mathrm{curl}\mathrm{curl}{\mathbf{u}}_{\epsilon }+\mathrm{div}({\mathbf{u}}_{\epsilon }\otimes \mathbf{\beta })+\sigma {\mathbf{u}}_{\epsilon }+\nabla p_{\epsilon }=$$

$$\mathrm{div}{\mathbf{u}}_{\epsilon }=$$

$${\mathbf{u}}_{\epsilon }\cdot \mathbf{n}=$$

$$\mathrm{curl}{\mathbf{u}}_{\epsilon }\times \mathbf{n}=$$

$$\epsilon (\mathrm{curl}{\mathbf{u}}_{\epsilon },\mathrm{curl}\mathbf{v})-({\mathbf{u}}_{\epsilon },\mathbf{\beta }\cdot \nabla \mathbf{v})+(\sigma {\mathbf{u}}_{\epsilon },\mathbf{v})-(p_{\epsilon },\mathrm{div}\mathbf{v})=$$

$$(\mathrm{div}{\mathbf{u}}_{\epsilon },q)=$$

$$\sqrt{\epsilon }\Vert \mathrm{curl}{\mathbf{u}}_{\epsilon }{\Vert }_{\Omega }+\Vert \sqrt{\sigma }{\mathbf{u}}_{\epsilon }{\Vert }_{\Omega }+\Vert p_{\epsilon }{\Vert }_{\Omega }\le C\Vert \mathbf{f}{\Vert }_{\Omega }.$$

Full text

1. Introduction

The use of H(div)-conforming finite element methods for the approximation of incompressible flow at high Reynolds number has been receiving increasing attention from the research community recently [14, 20, 25]. By construction such methods can satisfy the divergence-free condition exactly. The lack of $H^{1}$-conformity is handled using techniques drawing on ideas from discontinuous Galerkin methods [12], resulting in several possible different choices for the discretisation of the transport term and the viscous term. For the former one may either design an energy conserving method using central fluxes, or one may opt for a dissipative alternative in the form of upwind fluxes. The latter were shown in [14] to be more robust than the former, as is the case for discontinuous Galerkin (DG) methods. For DG-methods applied to scalar problems it is well known that thanks to the dissipative properties of the upwind flux one may prove an error estimate in the $L^{2}$-norm, of the form (see, e.g., [17])

[TABLE]

where $u$ is the exact solution, $u_{h}$ its DG-approximation, $\Omega\subset\mathbb{R}^{d}$, $d=2,3$. is the computational domain, $h$ the mesh parameter, and finally $k$ the polynomial degree of the approximation space. On special meshes one can in fact prove optimal estimates with rate $h^{k+1}$ for upwind DG methods applied to scalar problems [6, 23]. However, as it is shown in [21], the result (1.1) is sharp on general meshes.

Estimates of the type (1.1) are also the best that are known for either stabilised conforming finite element approximations, or fully DG methods, of laminar solutions of the Navier-Stokes’ equations in the high Reynolds number regime [5, 15], or the incompressible Euler equations [18, 4]. The robustness of the H(div)-conforming elements in the case of vanishing viscosity was shown in [19] for the case of the Brinkman problem, i.e. without the convection terms. Despite all the work quoted above, there seems to be no proof of an error estimate of the type (1.1) for finite element methods using H(div)-conforming elements applied to incompressible flow problems (see the discussion in [25, 20]).

The purpose of this work is to fill the gap mentioned in the last paragraph. That is, proving an estimate of the type (1.1) for finite element methods approximating a stationary linearised model of inviscid flow and using H(div)-conforming approximation spaces for the velocity approximation. Both the spaces designed by Raviart and Thomas [22] and by Brezzi, Douglas and Marini [3] enter the framework. As stabilising fluxes, these need to be either upwind, or, in case of central fluxes, an additional penalty term on the jump of the tangential component of the velocity needs to be added. In the particular case in which the velocity is approximated using the Raviart-Thomas space we also prove a convergence result for the pressure error, showing that the approximate pressure converges to the exact pressure in the $L^{2}$-norm also with the rate $O(h^{k+\frac{1}{2}})$. For the BDM space the rate $O(h^{k+\frac{1}{2}})$ is obtained for the projection of the error onto the pressure space, but since in this case the pressure space is of polynomial degree $k-1$, this is a superconvergence result.

1.1. Linear model problem

To keep the discussion as simple as possible we consider the following linear model problem.

Find a velocity $\bm{u}$ and a pressure $p$ satisfying

[TABLE]

We think of $\bm{u}$ and ${\bm{\beta}}$ as column vectors and we set $\bm{u}\otimes{\bm{\beta}}=\bm{u}{\bm{\beta}}^{t}$. We assume that ${\mathop{\mathrm{div}\,}}{\bm{\beta}}=0$ and that $\sigma\in L^{\infty}(\Omega)$ with $\sigma(\bm{x})\geq\sigma_{0}>0$ almost everywhere in $\Omega$. We assume that ${\bm{\beta}}\cdot\bm{n}=0$ on $\Gamma$. In spite of it being the natural candidate for a model problem for the development and analysis of numerical methods for inviscid flow this model does not seem to have been considered in the literature. Below we will first discuss the flow modelling leading to the system (1.2).

To obtain the stationary linear model problem (1.2) from the incompressible Euler equations, assume that a stationary solution to the latter ${\bm{\beta}}$, is subject to a smooth, exponentially growing perturbation of the right hand side of the momentum equation of the form:

[TABLE]

Writing the perturbed solution ${\bm{\beta}}+\tilde{\bm{u}}$ where $\tilde{\bm{u}}(x,t)$ is the perturbation resulting from the pertubation of the right hand side and neglecting quadratic terms in the perturbation $\tilde{\bm{u}}$, we may write the linearised momentum equation

[TABLE]

With the above choice of perturbation we may write the solution on the separated form

[TABLE]

Injecting this expression in (1.3) we arrive at the following stationary form for the space varying part of the perturbation

[TABLE]

To further simplify the model problem we finally drop the second term in the left hand side of (1.4). Since ${\mathop{\mathrm{div}\,}}({\bm{\beta}}\otimes\bm{u})=\bm{u}\cdot\nabla{\bm{\beta}}$, this is a non-essential term which can be absorbed in the reaction term under suitable assumptions on $\sigma$ and ${\bm{\beta}}$.

It is easy to construct solutions to the system (1.2). Examples of such solutions in the unit square are

1. (1)

x-independent solution.

Let ${\bm{\beta}}\cdot\bm{n}=0$ on $y=0$ and $y=1$ and ${\bm{\beta}}$ is defined to be periodic at $x=0$ and $x=1$. Then for any function $\varphi:\mathbb{R}\mapsto\mathbb{R}$, $\varphi\in[C^{1}(\mathbb{R})]^{2}$ a solution is given by:

[TABLE]

The associated pressure is $p=0$. 2. (2)

Stationary vortex sheet.

Let ${\bm{\beta}}\cdot\bm{n}=0$ on the boundaries of the square and define the streamfunction $\varphi(x,y):=\sin(n\pi x)\sin(n\pi y)$, corresponding to the vorticity $\omega:=\Delta\varphi=-2n^{2}\pi^{2}\sin(n\pi x)\sin(n\pi y)=-2n^{2}\pi^{2}\varphi$ with $n$ a positive integer. Then define:

[TABLE]

Since ${\bm{\beta}}\cdot\nabla\omega=-2n^{2}\pi^{2}(\partial_{y}\varphi(x,y)\partial_{x}\varphi(x,y)-\partial_{x}\varphi(x,y)\partial_{y}\varphi(x,y))=0$ we see that ${\bm{\beta}}$ is a solution to the two-dimensional stationary equations of inviscid flow. It is straightforward to verify that the velocity pressure formulation is satisfied for the pressure,

[TABLE]

In both examples (1) and (2) we achieve a problem on the form (1.2) by taking ${\bm{f}}=\sigma{\bm{\beta}}$ and the solution is then $\bm{u}={\bm{\beta}}$.

1.2. Outline of paper

We prove existence of solutions of the model problem (1.2) and uniqueness for $\sigma$ large enough, on smooth domains, in section 3. The H(div)-conforming upwind finite element methods are introduced and analysed in section 4. Finally in section 5 we illustrate the theory by computing approximations to the example (2) above.

2. Notation and preliminary results

The partial differential equation will be posed on an open polyhedral domain $\Omega\subseteq\mathbb{R}^{d},d=2,3$ with Lipschitz boundary $\Gamma$. For some of the theoretical results we will assume a smoother boundary. We adopt standard notation for Sobolev and Lebesgue spaces. In particular, for $D\subset\Omega$ we denote by $(\cdot,\cdot)_{D}$ the $L^{2}(D)$ inner product (without making a distinction between scalar and vector and tensor-valued functions). For $D=\Omega$ we drop the subindex in the above notation. The norm in $L^{2}(D)$ will be denoted by $\|\cdot\|_{D}$. By $W^{m,p}(D),m\geq 0,1\leq p\leq\infty$ we will denote the functions in $L^{p}(D)$, with distributional derivatives up to order $m$ belonging to $L^{p}(D)$, with norm (seminorm) $\|\cdot\|_{m,p,D}$ ($|\cdot|_{m,p,D}$). For $p=2$ we denote $H^{m}(D)=W^{m,2}(D)$, and the corresponding norm is denoted $\|\cdot\|_{m,D}$. As usual, $H^{m}_{0}(D)$ denotes the closure of $C^{\infty}_{0}(D)$ in the $\|\cdot\|_{m,D}$-norm. We also denote by $L^{2}_{0}(D)$ the space of $L^{2}(D)$ functions with zero mean value in $D$. All spaces for vector-valued functions will be denoted by boldface notation, e.g., $\bm{H}^{1}(D)=[H^{1}(D)]^{d}$, hence we denote by $\bm{H}({\mathop{\mathrm{div}\,}},D)$ the space of $\bm{L}^{2}(D)$ functions with distributional divergence in $\bm{L}^{2}(D)$, $\bm{H}_{0}({\mathop{\mathrm{div}\,}},D)=\{\bm{v}\in\bm{H}({\mathop{\mathrm{div}\,}},D):\bm{v}\cdot\bm{n}=0\;\textrm{on}\;\partial D\}$, and $\bm{H}({\mathop{\mathrm{curl}\,}},D)$ denotes the space of $\bm{L}^{2}(D)$ functions with distribution curl in $\bm{L}^{2}(D)$.

Below we will make use of the following preliminary result (for its proof, see, e.g., [11]).

Lemma 2.1**.**

There exists a constant $C>0$ such that for every $q\in L_{0}^{2}(\Omega)$ there exists $\bm{v}\in\bm{H}_{0}^{1}(\Omega)$ satisfying

[TABLE]

[TABLE]

Also in [11] the proof of the following result can be found.

Proposition 2.2**.**

The following bound holds

[TABLE]

If we assume that $\partial\Omega$ is $C^{1,1}$

[TABLE]

Finally, if $\Omega$ is a convex Lipschitz polyhedron [1], or a convex more regular domain, then

[TABLE]

Finally, for two $3\times 3$ matrices $A$ and $B$ with rows $A_{i}$ and $B_{i}$ ($i=1,2,3$) we define $C:=A\times B$ with $C_{1}=A_{2}\cdot B_{3}-A_{3}\cdot B_{2}$, $C_{2}=-(A_{1}\cdot B_{3}-A_{3}\cdot B_{1})$ $C_{3}=A_{1}\cdot B_{2}-A_{2}\cdot B_{1})$, and a simple calculation gives the following identity.

Lemma 2.3**.**

It holds

[TABLE]

3. Well-Posedness of the model problem

It appears that the linear inviscid model (1.2) has not been analysed mathematically. Hence, will here first study its well-posedness before proceeding with the finite element analysis. Transport problems have been studied by several authors (e.g. [13, 10, 8]). However, the incompressibility constraint seems to add new challenges to the analysis and we cannot apply the techniques of the above mentioned papers directly. The weak formulation of (1.2) is given by:

Find $\bm{u}\in\bm{H}_{0}({\mathop{\mathrm{div}\,}},\Omega)$ and $p\in L_{0}^{2}(\Omega)$ that satisfy

[TABLE]

3.1. Existence of weak solutions (3.1)

In order to prove existence of the problem (3.1) we will regularize it. Consider the following problem: Find a velocity $\bm{u}_{\varepsilon}$ and a pressure $p_{\varepsilon}$ satisfying

[TABLE]

The weak formulation of (3.2) is as follows: Find $(\bm{u}_{\varepsilon},p_{\varepsilon})\in\bm{H}_{0}^{1}(\Omega)\times L_{0}^{2}(\Omega)$ such that

[TABLE]

Lemma 3.1**.**

There exists a unique solution $\bm{u}_{\varepsilon}\in\bm{H}_{0}^{1}(\Omega)$ and $p_{\varepsilon}\in L_{0}^{2}(\Omega)$ to the problem (3.3). In addition, if ${\bm{\beta}}\in\bm{L}^{\infty}(\Omega)$, then the following bound holds

[TABLE]

where the constant $C$ depends on $\sigma$ and $\|{\bm{\beta}}\|_{\infty,\Omega}$, but not on negative powers of $\varepsilon$.

Proof.

Existence and uniquness of a solution of (3.3) follows from the Babuska-Brezzi theory [2]. Testing the equation with $\bm{u}_{\varepsilon}$ we get

[TABLE]

Therefore, we have the bound

[TABLE]

Moreover, using Lemma 2.1 and (3.3a) we have that

[TABLE]

and the proof is finished using (3.6). ∎

Theorem 3.2**.**

There exists a solution $\bm{u}\in\bm{L}^{2}(\Omega)$ and $p\in L^{2}(\Omega)$ to (3.1).

Proof.

Since $\{\bm{u}_{\varepsilon}\}$ and $\{p_{\varepsilon}\}$ are uniformly bounded in $\bm{H}_{0}({\mathop{\mathrm{div}\,}},\Omega)$ and $L^{2}_{0}(\Omega)$, respectively, there exists a subsequence such that $\bm{u}_{\varepsilon}\rightharpoonup\bm{u}$ and $p_{\varepsilon}\rightharpoonup p$ with $\bm{u}\in\bm{H}_{0}({\mathop{\mathrm{div}\,}},\Omega)$ and $p\in L_{0}^{2}(\Omega)$. Moreover, since ${\mathop{\mathrm{div}\,}}\bm{u}_{\varepsilon}=0$, for all $\phi\in H_{0}^{1}(\Omega)$ we have $(\bm{u},\nabla\phi)=\lim_{\varepsilon\rightarrow 0}(\bm{u}_{\varepsilon},\nabla\phi)=\lim_{\varepsilon\rightarrow 0}-({\mathop{\mathrm{div}\,}}\bm{u}_{\varepsilon},\phi)=0$, thus showing that ${\mathop{\mathrm{div}\,}}\bm{u}=0$ in $\Omega$. We then see that from (3.3a) and the fact that $\varepsilon\|\nabla\bm{u}_{\varepsilon}\|_{\Omega}\rightarrow 0$ as $\varepsilon\to 0$ that $\bm{u}$ and $p$ satisfy (3.1). ∎

3.2. Uniqueness of weak solutions

In general we cannot prove uniqueness of weak solutions (3.1). However, we will be to prove existence and uniqueness of solutions in the space $\bm{H}^{1}(\Omega)\cap\bm{H}_{0}({\mathop{\mathrm{div}\,}},\Omega)$ by making more stringent requirements on the coefficients and the boundary $\Gamma$. To achieve this goal, it is necessary to introduce a different regularised (as compared to (3.2)) problem to prove existence of smoother solutions to (3.1). The idea consists in considering the folllowing regularised Hodge-Oseen problem: Find a velocity $\bm{u}_{\varepsilon}$ and a pressure $p_{\varepsilon}$ satisfying

[TABLE]

The weak formulation of (3.7) reads as follows: Find $\bm{u}_{\varepsilon}\in\bm{V}:=\bm{H}^{1}(\Omega)\cap\bm{H}_{0}({\mathop{\mathrm{div}\,}},\Omega)$ and $p_{\varepsilon}\in L_{0}^{2}(\Omega)$ that satisfy

[TABLE]

Theorem 3.3**.**

Assume that $\bm{f}\in\bm{L}^{2}(\Omega)$ and that $\Gamma$ is $C^{1,1}$, or $\Omega$ is a convex Lipschitz polyhedron. Then, there exists a unique solution of (3.8). In addition, it satisfies

[TABLE]

Moreover, suppose that $\bm{f}\in\bm{H}^{1}(\Omega),{\bm{\beta}}\in\bm{W}^{1,\infty}(\Omega)$, $\sigma\in W^{1,\infty}(\Omega)$ and $\Gamma$ is $C^{3}$. If $\Omega$ is convex, let $\mathcal{C}=\|\nabla{\bm{\beta}}\|_{L^{\infty}(\Omega)}$, or otherwise $\mathcal{C}=K\|\nabla{\bm{\beta}}\|_{\infty,\Omega}$ where $K$ is from (2.2). Then, assuming $\sigma_{0}>\mathcal{C}$ we have

[TABLE]

where $C>0$ depends on $\sigma,{\bm{\beta}}$, and $K$, but not on negative powers of $\varepsilon$.

Proof.

The existence and uniqueness of this solution follows from the Babuska-Brezzi theory [2] by noting that as proven in [11], the norm in $\bm{H}^{1}(\Omega)$ is equivalent to the one in $\bm{H}({\mathop{\mathrm{curl}\,}},\Omega)\cap\bm{H}_{0}({\mathop{\mathrm{div}\,}},\Omega)$, thanks to the hypotheses on $\Gamma$. The bound (3.9) follows taking $\bm{v}=\bm{u}_{\varepsilon}$ in (3.8a), and the inf-sup conditions provides the stability for $p_{\varepsilon}$.

Next, whenever we suppose that $\Gamma$ is of class $C^{3}$ and $\bm{f}\in\bm{H}^{1}(\Omega)$, using the results in [26] (see Thereom 12 and Remark 16) we have the regularity $\bm{u}_{\varepsilon}\in\bm{H}^{3}(\Omega)$ and $p_{\varepsilon}\in H^{2}(\Omega)$. Noting that ${\mathop{\mathrm{curl}\,}}\bm{u}_{\varepsilon}\times\bm{n}=0$ on $\Gamma$ it follows that ${\mathop{\mathrm{curl}\,}}{\mathop{\mathrm{curl}\,}}(\bm{u}_{\varepsilon})\cdot\bm{n}=0$, so, $\tilde{\bm{v}}:={\mathop{\mathrm{curl}\,}}{\mathop{\mathrm{curl}\,}}(\bm{u}_{\varepsilon})\in\bm{H}^{1}(\Omega)\cap\bm{H}_{0}({\mathop{\mathrm{div}\,}},\Omega)$, and then it is a valid test function to be used in (3.8). Thus, taking $\tilde{\bm{v}}$ as test function in (3.8) and integrating by parts we obtain

[TABLE]

The second term in the left can be written as

[TABLE]

However, using Lemma 2.3 and the antisymmetry of the convective term

[TABLE]

and then

[TABLE]

Therefore, replacing the last identity in (3.10) we have

[TABLE]

Using the Cauchy Schwarz inequality, one of the inequalities (2.2) or (2.3), and the fact that ${\mathop{\mathrm{div}\,}}\bm{u}_{\varepsilon}=0$ we have

[TABLE]

Hence,

[TABLE]

and the proof follows dividing by $\|{\mathop{\mathrm{curl}\,}}\bm{u}_{\varepsilon}\|_{\Omega}$ and applying (3.9). ∎

Theorem 3.4**.**

Let us assume the hypotheses Theorem 3.3 . Then, there exists a unique solution of (3.1) such that $\bm{u}\in\bm{H}^{1}(\Omega)\cap\bm{H}_{0}({\mathop{\mathrm{div}\,}},\Omega)$ and $p\in H^{1}(\Omega)$.

Proof.

Let $(\bm{u}_{\varepsilon},p_{\varepsilon})$ be the solution of (3.8). Then, by Theorem 3.3 $\{(\bm{u}_{\varepsilon},p_{\varepsilon})\}$ is uniformly bounded in $\bm{H}^{1}(\Omega)\times L^{2}_{0}(\Omega)$. Hence, there exists a subsequence such that $\bm{u}_{\varepsilon}\rightharpoonup\bm{u}\in\bm{H}^{1}(\Omega)$ and $p_{\varepsilon}\rightharpoonup p\in L_{0}^{2}(\Omega)$ weakly. Moreover, since ${\mathop{\mathrm{div}\,}}\bm{u}_{\varepsilon}=0$ and $\bm{u}_{\varepsilon}\in\bm{H}_{0}({\mathop{\mathrm{div}\,}},\Omega)$, then ${\mathop{\mathrm{div}\,}}\bm{u}=0$ and $\bm{u}\in\bm{H}_{0}({\mathop{\mathrm{div}\,}},\Omega)$. This proves that $\bm{u}$ satisfies the second equation in (3.1) and the boundary conditions. Since $\bm{u}_{\varepsilon}\rightharpoonup\bm{u}\in\bm{H}^{1}(\Omega)$ weakly, then $\bm{u}_{\varepsilon}\to\bm{u}$ strongly in $\bm{L}^{2}(\Omega)$. In addition, since $\varepsilon\|{\mathop{\mathrm{curl}\,}}\bm{u}_{\varepsilon}\|_{\Omega}\to 0$ as $\varepsilon\to 0$, then using the weak convergence of $p_{\varepsilon}$ to $p$ in $L^{2}(\Omega)$ we can take the limit as $\varepsilon\to 0$ in (3.8) and conclude that $(\bm{u},p)$ also satisfies the first equation in (3.1). Finally, from the first equation in (3.1) we have $\nabla p=\bm{f}-\sigma\,\bm{u}-{\mathop{\mathrm{div}\,}}({\bm{\beta}}\otimes\bm{u})\in\bm{L}^{2}(\Omega)$, and then $p\in H^{1}(\Omega)$.

To prove uniqueness, assume that $\bm{f}=\bm{0}$. If we test with $\bm{u}\in\bm{H}^{1}(\Omega)\cap\bm{H}_{0}({\mathop{\mathrm{div}\,}},\Omega)$ we immediately get that $\|\sqrt{\sigma}\bm{u}\|_{\Omega}^{2}=0$ which gives that $\bm{u}=\bm{0}$. It easily follows that $p=0$. ∎

We finish this section by stating the following result that, in essence, casts the problem (3.1) as the limit of the Oseen problem (3.3).

Corollary 3.5**.**

Under the same hypotheses from Theorem 3.4 the solution $(\bm{u},p)$ of (3.1) is the limit of the solutions of the Oseen problem (3.3) in the following sense

[TABLE]

Proof.

The error $(\bm{u}-\bm{u}_{\varepsilon},p-p_{\varepsilon})$ satisfies the following error equation

[TABLE]

for all $\bm{v}\in\bm{H}^{1}(\Omega)\cap\bm{H}_{0}({\mathop{\mathrm{div}\,}},\Omega)$. Since $\bm{u}\in\bm{H}^{1}(\Omega)\cap\bm{H}_{0}({\mathop{\mathrm{div}\,}},\Omega)$, $\hat{\bm{v}}:=\bm{u}-\bm{u}_{\varepsilon}$ is a valid test function for (3.12). So, using $\hat{\bm{v}}$ in (3.12), the fact that both $\bm{u}$ and $\bm{u}_{\varepsilon}$ are divergence-free, the Cauchy-Schwarz inequality, and (3.4) we get

[TABLE]

as $\varepsilon\to 0$, which proves the convergence of $\bm{u}_{\varepsilon}$ to $\bm{u}$ in $\bm{L}^{2}(\Omega)$. The convergence of $\bm{u}_{\varepsilon}$ to $\bm{u}$ in $\bm{H}_{0}({\mathop{\mathrm{div}\,}},\Omega)$ follows from the fact that both $\bm{u}_{\varepsilon}$ and $\bm{u}$ are divergence-free.

To prove the convergence of the pressure, using Lemma 2.1 there exists $\bm{w}\in\bm{H}^{1}_{0}(\Omega)$ such that $|\bm{w}|_{1,\Omega}\leq C\|p-p_{\varepsilon}\|_{\Omega}$ and ${\mathop{\mathrm{div}\,}}\bm{w}=p-p_{\varepsilon}$. Then, using (3.12), the Cauchy-Schwarz inequality, and the convergence of $\bm{u}_{\varepsilon}$ to $\bm{u}$,

[TABLE]

and the proof follows by dividing by $\|p-p_{\varepsilon}\|_{\Omega}$ and noticing that, thanks to (3.4) the term within parentheses tends to zero as $\varepsilon\to 0$. ∎

4. Upwind H(div) method

4.1. Preliminaries

We denote by $\{\mathcal{T}_{h}\}_{h>0}$ a family of shape-regular simplicial triangulations of $\Omega$. The elements of $\mathcal{T}_{h}$ are denoted by $T$, with diameter $h_{T}$, and $h:=\max\{h_{T}:T\in\mathcal{T}_{h}\}$. The set of its facets (edges for $d=2$, faces for $d=3$) is denoted by $\mathcal{E}_{h}$. To cater for the nonconforming character of the approximation we also introduce the following broken versions of the scalar product

[TABLE]

In addition, we introduce the broken space $H(\mathcal{T}_{h})$, of functions in $L^{2}(\Omega)$ whose restriction to every $T\in\mathcal{T}_{h}$ belongs to $H(T)$.

Let $T\in\mathcal{T}_{h}$ and let $\bm{x}\in\partial T$ then we define

[TABLE]

and

[TABLE]

For $F\in\mathcal{E}_{h}$ and $F=\partial T_{1}\cap\partial T_{2}$ for $T_{1},T_{2},\in\mathcal{T}_{h}$ we define the jumps

[TABLE]

and for $F\in\mathcal{E}_{h}$ and $F\subset\Gamma$ we define

[TABLE]

We then define the semi-norm on the jumps of the solution over element boundaries to be

[TABLE]

With these definitions we can state the following important identity [12, Lemma 6.1]

Proposition 4.1**.**

For all $\bm{v}\in\bm{H}^{1}(\mathcal{T}_{h})$, the following holds

[TABLE]

Let us define the Raviart-Thomas [22] and BDM spaces [3]. The space of polynomials of degree at most $k$ defined in $T$ is denoted by $\EuScript{P}_{k}(T)$, and we denote $\bm{\EuScript{P}}_{k}(T)=[\EuScript{P}_{k}(T)]^{d}$. For every $T\in\mathcal{T}_{h}$, let $\text{RT}_{k}(T)=\bm{\EuScript{P}}_{k}(T)+(\EuScript{P}_{k}(T)\setminus\EuScript{P}_{k-1}(T))\bm{x}$. We define, for $k\geq 0$, the spaces

[TABLE]

A well-known property linking these two spaces is stated now (for a proof see [7, Lemma 4.3]).

Lemma 4.2**.**

Let $\bm{v}\in\bm{V}_{h,k}^{\text{RT}}$ with ${\mathop{\mathrm{div}\,}}\bm{v}=0$ on $\Omega$ then $\bm{v}\in\bm{V}_{h,k}^{\text{BDM}}$.

We next introduce the standard $L^{2}$-projection on polynomials on an element $T$, $P_{k}^{T}:L^{2}(T)\rightarrow\EuScript{P}_{k}(T)$. Its global equivalent will be denoted $P_{k}:L^{2}(\Omega)\rightarrow M_{h,k}$. We recall the standard estimates for the $L^{2}$-projection (see, e.g., [9])

[TABLE]

The Raviart-Thomas interpolation operator will be used in the sequel. It is defined as follows: $\bm{\Pi}:\bm{H}^{1}(\Omega)\cap\bm{H}_{0}(\text{div};\Omega)\rightarrow\bm{V}_{h,k}^{\text{RT}}$ where $\bm{\Pi}\bm{v}$ is the only function of $\bm{V}_{h,k}^{RT}$ satisfying

[TABLE]

This operator satisfies the following classical properties (see, e.g., [2]).

Lemma 4.3**.**

Let $k\geq 0$. The mapping $\Pi$ satisfies the following commutative property

[TABLE]

Let $\bm{v}\in\bm{H}^{k+1}(\Omega)$ then we have

[TABLE]

We end this section recalling the following classical inverse and local trace inequalities that hold for every $T\in\mathcal{T}_{h}$

[TABLE]

4.2. The finite element method and the error estimates for the velocity

Throughout, the velocity and pressure will be approximated using the spaces $\bm{V}_{h}$ and $M_{h}$, respectively. In this work we will consider the following choices:

[TABLE]

or

[TABLE]

The numerical method analysed here reads: Find $\bm{u}\in\bm{V}_{h}$ and $p_{h}\in M_{h}$ such that

[TABLE]

Thanks to the inf-sup stability of the pair $\bm{V}_{h}\times M_{h}$ (see [2]), and Proposition 4.1, problem (4.9) has a unique solution. Moreover, the method (4.9) is consistent; in fact, for $(\bm{u},p)\in\bm{H}^{1}(\Omega)\times L^{2}_{0}(\Omega)$ solving (1.2) we have

[TABLE]

A consequence of Lemma 4.2 is that the finite element method (4.9) produces the same velocity approximation for $\bm{u}_{h}\in\bm{V}_{h,k}^{\text{RT}}$ and $\bm{u}_{h}\in\bm{V}_{h,k}^{\text{BDM}}$. We show that in the following proposition.

Proposition 4.4**.**

Let $(\bm{u}_{h},p_{h})$ be the solution of (4.9) for the spaces $\bm{V}_{h}\times M_{h}=\bm{V}_{h,k}^{\text{RT}}\times M_{h,k}$ and $(\tilde{\bm{u}}_{h},\tilde{p}_{h})$ the solution of (4.9) for the spaces $\bm{V}_{h}\times M_{h}=\bm{V}_{h,k}^{\text{BDM}}\times M_{h,k-1}$. Then $\bm{u}_{h}=\tilde{\bm{u}}_{h}$.

Proof.

Let $\bm{e}_{h}:=\tilde{\bm{u}}_{h}-\bm{u}_{h}$, $\eta_{h}=\tilde{p}_{h}-p_{h}$ then using (4.9) we see that

[TABLE]

Since ${\mathop{\mathrm{div}\,}}\bm{e}_{h}=0$ by Lemma 4.2there holds $\bm{e}_{h}\in\bm{V}_{h,k}^{\text{BDM}}$, which is a valid test function. Taking $\bm{v}_{h}=\bm{e}_{h}$ in (4.11) and applying Proposition 4.1 we obtain

[TABLE]

which proves the claim. ∎

We can now derive an error estimate for the velocity. We let $\bm{e}_{h}=\bm{\Pi}\bm{u}-\bm{u}_{h}$ and start by noticing that

[TABLE]

Hence, by Lemma 4.2 we have $\bm{e}_{h}\in\bm{V}_{h,k}^{\text{BDM}}$ and in particular

[TABLE]

Theorem 4.5**.**

Let $\bm{u}\in[H^{1}(\Omega)]^{d}$ solve (1.2) and let $\bm{u}_{h}\in\bm{V}_{h}$ solve (4.9). Then, the following error estimate holds

[TABLE]

where the constant $C$ does not depend on $h$, or any physical parameter of the equation.

Proof.

Using (4.9), (4.10), (4.12), and (4.1) we get

[TABLE]

Hence, we have

[TABLE]

We bound each term separately. Using (4.13), the definition of $\bm{\Pi}$ (4.4)-(4.5), (4.3), and (4.7), we have

[TABLE]

Using the contributions from neighbouring elements on the face to express the discrete error on the faces in terms of jumps, the normal continuity of $\bm{u}$ and $\bm{\Pi}\bm{u}$, and using the local trace inequality (4.8) it is easy to show that

[TABLE]

Finally,

[TABLE]

Therefore, inserting (4.15)-(4.17) into (4.14) we arrive at

[TABLE]

The result follows after applying the triangle inequality. ∎

The following result appears as a corollary of the last theorem and Lemma 4.3.

Corollary 4.6**.**

Let $\bm{u}\in[H^{k+1}(\Omega)]^{d}$ solve (1.2) and let $\bm{u}_{h}\in\bm{V}_{h}$ solve (4.9). Then, the following error estimate holds

[TABLE]

Remark 4.7*.*

The arguments of Theorem 4.5 and Corollary 4.6 may be used to improve the order obtained Theorem 2.2 of [14] to $O(h^{k+\frac{1}{2}})$, if an upwind flux is used. Following the ideas above, use integration by parts in the first term of $I_{1}$ in the equation after (2.12). Then add and subtract the exact solution to the approximate solution in term $I_{3}$ and recombine terms, so that one may use continuity on the norm augmented with $L^{2}$-control on the faces the jumps of the approximate velocity.

4.3. $L^{2}$-error estimates for the pressure approximation

Since the pressure space is of polynomial degree $k$ for the method using the $RT$ space for velocity approximation and $k-1$ for the method using the $BDM$ space, the optimal order that can be obtained for the error of the pressure approximation in the $L^{2}$-norm is $O(h^{k+1})$ and $O(h^{k})$, respectively. Here we will prove the following orders for the pressure error :

1. (1)

in the first case (RT), $O(h^{k+\frac{1}{2}})$; this is, the same suboptimality of $O(h^{\frac{1}{2}})$ as for the velocity approximation. 2. (2)

in the second case (BDM) we get the optimal convergence $O(h^{k})$; considering that the pressure space is of degree $k-1$. For the discrete error, i.e. the projection of the error on the space $M_{h}$, we get an $O(h^{k+\frac{1}{2}})$ estimate, this is a superconvergence of $O(h^{\frac{1}{2}})$ compared with the approximation property of the space of constant functions.

Theorem 4.8**.**

Let $(\bm{u},p)\in\bm{H}^{1}(\Omega)\times L^{2}_{0}(\Omega)$ solve (1.2) and let $(\bm{u}_{h},p_{h})\in\bm{V}_{h}\times M_{h}$ solve (4.9). Let $\ell$ denote the polynomial degree of the space $M_{h}$. Then, the following error estimate holds

[TABLE]

Proof.

Using the surjectivity of the divergence operator as a mapping from $\bm{H}^{1}_{0}(\Omega)$ to $L^{2}_{0}(\Omega)$ there exists $\bm{v}_{p}\in\bm{H}^{1}_{0}(\Omega)$ such that ${\mathop{\mathrm{div}\,}}\bm{v}_{p}=P_{\ell}p-p_{h}$ and

[TABLE]

It follows from (4.18) and (4.6) that

[TABLE]

If $\bm{V}_{h}\equiv\bm{V}_{h,k}^{\text{RT}}$ then choose $\tilde{\bm{\Pi}}\bm{v}_{p}\in\bm{V}_{h,k}^{\text{RT}}$ and if $\bm{V}_{h}\equiv\bm{V}_{h,k}^{\text{BDM}}$ choose $\tilde{\bm{\Pi}}\bm{v}_{p}\in\bm{V}_{h,{k-1}}^{\text{RT}}\subset\bm{V}_{h,k}^{\text{BDM}}$. Using (4.9) and (4.10) we find that

[TABLE]

Applying the Cauchy-Schwarz inequality and the stability of the RT interpolant and of $\bm{v}_{p}$ we have

[TABLE]

For the remaining term observe that, by the definition of $\langle\cdot,\cdot\rangle_{h}$, the fact that ${\bm{\beta}}\cdot\bm{n}$ changes sign on neighbouring elements and that $(\bm{u}-\widehat{\bm{u}_{h}})$ is single valued on the faces of the triangulation,

[TABLE]

The right hand side of this equality is bounded using the Cauchy-Schwarz inequality, the trace inequality (4.8) and the interpolation properties of the RT-interpolant of Lemma 4.3 as follows

[TABLE]

where in the last step we added and subtracted $\bm{\Pi}\bm{u}$, used the triangle inequality and the inverse inequality (4.7). We conclude by using (4.18) . ∎

The following result is an immediate consequence of Theorem 4.8 and Corollary 4.6 and the approximation properties of the $L^{2}$-projection,

Corollary 4.9**.**

Assume that $\bm{V}_{h}=\bm{V}_{h,k}^{\text{RT}}$ and $M_{h}=M_{h,k}$. Then, there exists $\tilde{C}_{{\bm{\beta}},\sigma}>0$ that depends only on the constants in the bounds of Theorems 4.8 and Corollary 4.6 such that

[TABLE]

For the case in which $\bm{V}_{h}=\bm{V}_{h,k}^{\text{BDM}}$ and $M_{h}=M_{h,k-1}$, the following error estimate holds

[TABLE]

and

[TABLE]

where $\hat{C}_{{\bm{\beta}},\sigma}$ depends on the constants in the bounds of Theorems 4.8 and Corollary 4.6.

5. A numerical example

Here we will show some illustrations of the theory developed above using the analytical solution of example (2) in section 1.1. For ample qualitative numerical evidence of the performance of this type of method on physically relevant problems we refer to the references [25, 24].

We consider the domain $\Omega=(0,1)\times(0,1)$ and the solution (1.5)-(1.6) of example (2). We used the package FreeFEM++ [16] to implement the formulation (4.9) with either the BDM(1) element and piecewise constant pressures or the RT(1) element with piecewise affine, discontinuous, pressures. The linear systems were solved using UMFPACK and the meshes were of Union Jack type. In Tables 1-2 we report the errors of velocities and pressures in the (relative) $L^{2}$-norm. We also report the CPU time. We see that the velocity approximations have identical errors in the two cases as predicted by Proposition 4.4, whereas as expected the BDM(1) approximation has poorer convergence of the pressure. The RT(1) computation however is more costly by almost a factor three.

In Table 3 we report the variation of the error on a fixed mesh with $h=1/40$ and $\sigma=100$. The variable $n$, controlling the number of vortices, and hence influencing both $\|{\bm{\beta}}\|_{W^{1,\infty}(\Omega)}$ and $\|\bm{u}\|_{H^{2}(\Omega)}$ is taken in the set $n\in\{1,2,4,8\}$. We observe (approximately) linear growth in both velocity and pressures, except for the pressure for the method using the RT element, where the growth is stronger. For the highest value $n=8$, all errors are above $15\%$ on this mesh. In Table 4 we vary the coefficient $\sigma$ and see that also here the error growth for decreasing $\sigma$ is by and large linear for the velocities, as predicted by theory (Corollary 4.6) and the RT pressure (Corollary 4.9). The BDM pressure on the other hand is very robust with respect to variations in $\sigma$, but much larger than the RT-pressure. It starts increasing only for the smallest value of the parameter, when the pressure errors of the two approximation spaces are comparable. It follows that for small values of $\sigma$ the pressure approximation is of similar quality for the BDM and RT methods.

Acknowledgments

The work of Gabriel R. Barrenechea has been funded by the Leverhulme Trust through the Research Fellowship No. RF-2019-510. Erik Burman was partially supported by the grant: EP/P01576X/1. Johnny Guzman was partially supported by the grant: NSF, DMS # 1620100.