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?
  • $\begingroup$ 1) Multigrid 2) BFBt approximate commutators (see papers of Elman) $\endgroup$
    – Nick Alger
    Commented Apr 25, 2017 at 15:59
  • $\begingroup$ Typically either multigrid or possibly some sort of ILU for the A block. And you can approximate your scour complement with the plans matrix on the pressure space. Although, I would have a look into block preconditioners, they typically work much better. See 'Finite Elements and Fast iterative Solvers' by Elman, Silverster, Wathen. $\endgroup$
    – KyleW
    Commented Apr 28, 2017 at 7:31

2 Answers 2


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).

  • $\begingroup$ Do you know a C++-library which solves $(2)$ using such approaches? $\endgroup$
    – 0xbadf00d
    Commented Apr 25, 2017 at 21:07
  • $\begingroup$ Check out deal.II (link in the post); it may be useful to watch several lectures on it (#38 covers block preconditioners): math.colostate.edu/~bangerth/videos.html $\endgroup$
    – 56th
    Commented Apr 25, 2017 at 21:31
  • $\begingroup$ Also chapter 5 by Mardal et al. (diffpack.com/?f=book/overview) discusses how block preconditioners implemented in Dolfin. They discuss BD preconditioner for the Stokes problem there, but I believe it is possible extend this easily. $\endgroup$
    – 56th
    Commented Apr 25, 2017 at 21:57
  • $\begingroup$ In the paper of Rehman et al., are they suggesting to use $P$ (defined in $(11)$; with $F=A$) as a left preconditioner for $(2)$, i.e. solving $$P^{-1}\left(\begin{matrix}A&B^T\\B&0\end{matrix}\right)\left(\begin{matrix}u\\p\end{matrix}\right)=P^{-1}\left(\begin{matrix}f\\0\end{matrix}\right)\tag4\;,$$ or are they suggesting to replace $(2)$ by $$P\left(\begin{matrix}u\\p\end{matrix}\right)=\left(\begin{matrix}f\\0\end{matrix}\right)\tag5\;?$$ $\endgroup$
    – 0xbadf00d
    Commented Apr 25, 2017 at 22:00
  • $\begingroup$ They suggest to solve $(4)$. $\endgroup$
    – 56th
    Commented Apr 25, 2017 at 22:07

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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.