Preconditioner for the GMRES method in the Uzawa algorithm

Question

I'm trying to solve

\begin{equation}\left\{ \begin{split} \frac{\partial u}{\partial t}+(u\cdot\nabla)u-\nu\Delta u+\frac1\rho\nabla p&=f\;\;\;\text{in }\Lambda\\ u&=0\;\;\;\text{on }\partial\Lambda\\ \nabla\cdot u&=0\;\;\;\text{in }\Lambda \end{split}\right.\tag1 \end{equation}

on $\Lambda=(0,a)\times(0,b)$. Let

\begin{equation} \begin{split}\mathfrak a(u,v)&:=\sum_{i=1}^d\langle\nabla u_i,\nabla v_i\rangle_{L^2}\\ \mathfrak b(u,q)&:=-\langle\nabla\cdot u,q\rangle_{L^2}\\ \mathfrak c(u,v,w)&:=\langle(u\cdot\nabla)v,\cdot w\rangle_{L^2} \end{split} \end{equation}

for $u,v\in V:=H_0^1(\Lambda,\mathbb R^2)$ and $q\in\Pi:=\left\{p\in L^2(\Lambda):\int p=0\right\}$. Now, let

\begin{equation} \begin{split}\tilde{\mathfrak a}(u,v)&:=\langle u,v\rangle_{L^2}+\Delta t\nu\mathfrak a(u,v)+\Delta t\mathfrak c(u^0,u,v)\\ \tilde{\mathfrak b}(v,q)&:=\frac{\Delta t}\rho\mathfrak b(v,q) \end{split} \end{equation}

for $u,v\in V$ and $q\in\Pi$. $u^0$ is the approximate solution of the previous time step and $\Delta t$ is the elappsed time. I've chosen a Taylor-Hood pair for the finite element discretization. Now, I'm left with a system $$\left(\begin{matrix}A&B^T\\B&0\end{matrix}\right)\left(\begin{matrix}u\\p\end{matrix}\right)=\left(\begin{matrix}f\\0\end{matrix}\right)\tag2$$ where

\begin{equation} \begin{split} a_{ij}&=\tilde{\mathfrak a}(\phi_i,\phi_l)\\ b_{ij}&=\tilde{\mathfrak b}(\phi_j,\psi_i) \end{split} \end{equation}

and $(\phi_i)$ and $(\psi_i)$ are bases of the finite dimensional subspaces of $V$ and $\Pi$, respectively.

My idea was to solve

\begin{equation} \left\{\begin{split} BA^{-1}B^Tp&=BA^{-1}f\\ Au&=f-B^Tp\;. \end{split}\right.\tag3 \end{equation}

However, I've got several problems with $(3)$:

1. I'm trying to solve the linear equations using the GMRES method. As a first step, the right-hand side of the first equation, $BA^{-1}f$, has to be computed. I've assembled $A$ and $B$; hence I'm trying to find a solution $Ax=f$ (using the GMRES method) and then compute $Bx$. However, the GMRES method for $Ax=f$ converges extremely slowly (most probably due to the undesirable eigenvalue distribution of $A$). Can we accelerate the convergence by a preconditioner? If so, which one should I use?
2. The same question applies to the invocation of the GMRES method for $BA^{-1}B^Tp=BA^{-1}f$. Which preconditioner can I use here?

Answer 1

Please check this paper by Benzi et al. They address this issue and give corresponding references on p. 45.

Shortcut: for the Stokes problem $A = \text{diag}(A_{11},A_{11},\dots,A_{11})$ is just a collection of discrete Laplace operators, so it is natural to approximate their inverses using multigrid. However, things get much more complicated for Oseen type problems as in your case (especially if convection is high, $\nu \ll 1$); Uzawa converges rather slowly in this case. They give references to papers which propose preconditioning techniques for Uzawa.

Beyond Uzawa: it may be useful to look at different solving techniques, e.g. block preconditioners such as PCD preconditioner by Kay et al. or AL preconditioner by Benzi and Olshanskii. A nice overview is given by Rehman et al. A recent step of deal.II tutorial implements AL approach (yet they use a direct solver for $A$ for simplicity of implementation).

Answer 2

As already pointed out in the other answer, deal.II has a number of tutorial programs that deal with exactly this problem. In particular, step-22 uses essentially the approach you want to use, but concludes that this is among the least efficient ways to do things. In the appendix it points out a better way to do it, and provides a bit of code as well. It also has links to other tutorials such as step-57 that describe better and more efficient approaches.