Mixed Finite Element Method for the Stokes System—Some Implementation Details

Question

I am currently working on my bachelor’s diploma. The research concerns mixed finite element method for the 2D Stokes system

$$ - \Delta \boldsymbol u + \nabla p = \boldsymbol f, \quad \boldsymbol x \in \Omega \subset \mathbb R^2, \\ \nabla \cdot \boldsymbol u = 0, \quad \boldsymbol x \in \Omega, $$

with homogeneous Dirichlet (no–slip) BCs

$$ \boldsymbol u = \boldsymbol 0, \quad \boldsymbol x \in \partial \Omega. $$

Obtaining weak form and choosing appropriate spaces $\mathbb U_h \times \mathbb P_h \ni (\boldsymbol u_h, \, p_h)$ for the velocity field and pressure distribution, one reduces the problem to a saddle–point linear system

$$ \begin{bmatrix} \mathbf A & \mathbf B^T \\ \mathbf B & \mathbf 0 \end{bmatrix} \begin{bmatrix} \boldsymbol \xi \\ \boldsymbol \psi \end{bmatrix} = \begin{bmatrix} \boldsymbol b \\ \boldsymbol 0 \end{bmatrix}, $$

from which one then gets FE–solutions (using natural isomorphism $\mathcal P$ between vectors and FE–interpolants) for pressure $p = \mathcal P \, \boldsymbol \psi $ and ditto for velocity components.

I focus on the analysis and implementation of (geometric) multigrid solving techniques for this system.

I have an experience [Russian text] in solving elliptic and time–dependent parabolic and hyperbolic PDEs, yet I’ve never dealt with saddle–point problems before. So I want to clarify several aspects.

1. It is well–known that one should carefully choose FE–spaces in order to obtain a well–posed problem. I focus on two LBB–stable finite element pairs:

   i. $\Delta \, \boldsymbol P^1 \, CR — \Delta \, P^0 \, L$ finite element pair (non–conform space for the velocity components),

   ii. $\Delta \, \boldsymbol P^2 \, L — \Delta \, P^1 \, L$ finite element pair (also known as Taylor–Hood).

   Let me clarify these notations. I follow Ciarlet’s definition of FE here: $“\Delta”$ defines geometry (a triangle in my case), $“P^i”$ defines a polynomial space for shape functions, and the last component defines DOFs (which one uses to obtain shape functions and their cross–element behavior)—$“L”$ for “Lagrange” and $“CR”$ for “Crouzeix–Raviart.” I have some visualizations for (i) here.

   As far as I concerned, inf–sup stability of these elements should guarantee existence of the unique solution of the above linear system. However, since the momentum equation involves only pressure gradient, it is determined up to a constant; so one usually requires $p$ to have zero mean value: $\int_\Omega p \, d \boldsymbol x = 0$.

   The question is, should I enforce this constraint explicitly in the system? Authors of [1, p. 309] suggest the following modification (they use (i)^st FE pair):

   $$ \begin{bmatrix} \mathbf A & \mathbf B^T & \boldsymbol 0 \\ \mathbf B & \mathbf 0 & \boldsymbol a \\ \boldsymbol 0 & \boldsymbol a^T & 0 \end{bmatrix} \begin{bmatrix} \boldsymbol \xi \\ \boldsymbol \psi \\ \mu \end{bmatrix} = \begin{bmatrix} \boldsymbol b \\ \boldsymbol 0 \\ 0 \end{bmatrix}. $$Other publications I’ve seen so far do not mention such modifications. Moreover, since I am interested in implementation of geometric multigrid, I have to assemble prolongation / restriction operators (sparse matrices); this is pretty straightforward for velocity components and pressure vectors (however, non–conform case will require some special handling), yet it is not for $\mu$.

   UPD: I noticed this answer. It seems pretty reasonable that, since the “discrete version” of $\nabla p$ is $\mathbf B^T \boldsymbol \psi$, it should imply that operator $\mathbf B^T$ has a kernel consisting of the identity vector (I should check this by direct computations (for not–too–big matrices, of course) when I’m done with assembling routine).

   This also implies that the final system actually has an infinite number of solutions (doesn’t it?). However, this should not cause any problems for a Krylov solver to converge, and one can enforce the desired constraints on the pressure at the post–processing step.

   But now I got a bit confused. Does LBB–stability of a finite element pair really imply existence and uniqueness? Isn’t it really about existence?

2. According, for one, to [2, p. 99] it is reasonable to use multigrid as a preconditioner for a Krylov solver for div–grad elliptic problems. For the Stokes problem, is it reasonable to use multigrid (for one, with the Vanka smoother) as a preconditioner for MINRES? MINRES seems reasonable since the system is symmetric and sign–indefinite. I’ve seen articles suggesting using multigrid as a stand–alone solver.

3. I also want to compare my future implementation with “black–box” Krylov solvers (for one, in terms of residuals’ history). And I do not want to implement them all.

   I’ve already implemented data structures for CSC and symmetric CSlC sparse matrices (these data structures are used for $B$ and $A$ blocks of the system, respectively). I’ve also implemented import / export in Harwell–Boeing format for these data structures (HB–format is pretty popular and supported, for example, by Mathematica).

   So I am wondering which tools I can use. It would be nice if the solver routine provided residuals’ history (not only the solution vector). Any suggestions?

References:

1. Mats G. Larson, Fredrik Bengzon
   The Finite Element Method: Theory, Implementation, and Applications
   2013

2. Maxim A. Olshanskii, Eugene E. Tyrtyshnikov
   Iterative Methods for Linear Systems
   2014

Answer 1

The topic of linear solvers for the Stokes equation is pretty well discussed in the literature, and the common consensus is that (i) using GMRES as the outer solver, (ii) the Silvester-Wathen approach to forming a preconditioner for the outer solver, (iii) using multigrid as the inner solver or preconditioner for the elliptic block is the way to go. This is not to say that there are not other methods that can also be used, but this is what people consider state of the art.

My own take at implementing this is given in https://www.dealii.org/developer/doxygen/deal.II/step_22.html . The solver/preconditioner mentioned above is discussed in the "Possibilities for extensions section".

As for implementation: You are wasting your time if you try to implement everything yourself. Use PETSc or Trilinos to represent your matrices, vectors, and linear solvers. Use any of the finite element libraries that are out there to represent the mesh and everything finite element related. For example, you can write a simple Stokes solver with not very much more than 100 or 150 lines of code in deal.II (a project I am heavily involved with, see here).

Answer 2

Thanks a lot for the answers. I would like to summarize what have been said before and add some more details in this post. I hope it will help if anyone faces similar problems.

(1)

First note: indeed, inf–sup stability guarantees existence and uniqueness. Yet I missed a crucial point: it also guarantees optimal error estimates (thanks @knl). In particular, if your FE–pair is not stable, you may get nonphysical solutions, e.g. spurious modes of pressure (nonphysical oscillations) or locking phenomenon. Nice explanation and examples are given in [3, p. 62].

Second note: if we work w/ essential + natural BCs, we may choose spaces $\mathbb H^1(\Omega) \times \mathbb L^2(\Omega)$ and and we are fine. If we work w/ essential BCs only (as in my post), we may choose spaces $\mathbb H^1(\Omega) \times \mathbb L^2_0(\Omega)$ and we are fine again, at least mathematically. However, working w/ a discrete analogue of $\mathbb L^2_0(\Omega)$ is somewhat inconvenient in practice.

One way is to use approach from [1, p. 309] which I mentioned in my question. Another way is to fix one of the pressure DOFs (just like we usually handle essential BCs). A very good (and fun) explanation is given in [4, p. 293]:

As for my experience, I tried to solve the Stokes problem w/ essential BCs without any constraints on the pressure space. I used equivalent formulation where block (2, 1) equals $-\mathbf B$ instead of $\mathbf B$ (so the resulting matrix is not symmetric but sign–definite). I started with a random guess and used BiCGStab and PBiCGStab w/ block diagonal preconditioner. Convergence history in terms of number of iterations and working time were similar to the case w/ essential + natural BCs.

(2)

It is much more efficient and smart to use MG–iterations to approximate actions of inverses of diagonal blocks of $\mathbf A = \text{diag}(\mathbf A_{11}, \mathbf A_{22})$ since they are elliptic. This is a part in building physics–based / block preconditioners which are considered to be the most robust choice for “complex” problems if one wants to solve the system in $O(n)$ time.

Wolfgang Bangerth mentioned this approach. It was extremely helpful for me to watch his lectures on preconditioners (especially this one which covers block preconditioners), so I marked his answer.

I decided not to implement Vanka’s approach—I implemented block diagonal preconditioner for the problem in my question and AL–preconditioner for the Oseen problem.

References:

3. Volker John
   Finite Element Methods for Incompressible Flow Problems
   2016
4. Daniele Boffi, Franco Brezzi, Michel Fortin
   Mixed Finite Element
   Methods and Applications
   2013

Answer 3

I happen to work on the same mixed finite element method as well, but in a different context. (So I don't know if my answer can fully apply to your problem)

In my work I'm enforcing the compability constraint directly. I've seen the extended system of equations you have posted numerous times. In solid mechanics, it will be used as a kind of stabilization, but is not really needed, at least to my experience. A big downfall (of mixed methods in general) which might be able to overcome with the extension of the system of equations is that the mixed stiffnessmatrix (or however you want to call it...) is ill-conditioned. It resulted in a poor convergence in my case.