# Impose the compatibility conditions for mixed finite elements method in Stokes equation

$\newcommand{\v}[1]{\boldsymbol{#1}}$ Suppose we have following Stokes flow model equation:

\tag{1} \left\{ \begin{aligned} -\mathrm{div}(\nu \nabla \v{u}) + \nabla p &= \v{f} \\ \mathrm{div} \v{u} &= 0 \end{aligned} \right. where the viscosity $\nu(x)$ is a function, for the standard mixed finite element, say we use the stable pair: Crouzeix-Raviart space $\v{V}_h$ for the velocity $\v{u}$ and element-wise constant space $S_h$ for the pressure $p$, we have the following variational form:

$$\mathcal{L}([\v{u},p],[\v{v},q]) = \int_{\Omega} \nu \nabla\v{u}:\nabla \v{v} -\int_{\Omega} q\mathrm{div} \v{u} -\int_{\Omega} p\mathrm{div} \v{v} =\int_{\Omega} \v{f}\cdot \v{v} \quad \forall \v{v}\times q\in \v{V}_h\times S_h$$

and we know that since the Lagrange multiplier $p$ can be determined up to a constant, the finally assembled matrix should have nullspace $1$, to circumvent this we could enforce the pressure $p$ on some certain element be zero, so that we don't have to solve a singular system.

So here is my question 1:

• (Q1) Is there other way than enforcing $p=0$ on some element to eliminate the kernel for standard mixed finite element? or say, any solver out there that be able to solve the singular system to get a compatible solution?(or some references are welcome)

And about the compatibility, for (1) it should be $$\int_{\Omega} \nu^{-1} p = 0$$ and the nice little trick is to compute $\tilde{p}$ be the $p$ we got from the solution of the linear system subtracted by its weighted average: $$\tag{2} \tilde{p} = p - \frac{\nu}{|\Omega|}\int_{\Omega} \nu^{-1} p$$

However, recently I have just implemented a stabilized $P_1-P_0$ mixed finite element for Stokes equation by Bochev, Dohrmann,and Gunzberger, in which they added a stabilized term to the variational formulation (1): $$\widetilde{\mathcal{L}}([\v{u},p],[\v{v},q]) =\mathcal{L}([\v{u},p],[\v{v},q]) -\int_{\Omega} (p - \Pi_1 p)(q -\Pi_1 q) =\int_{\Omega} \v{f}\cdot \v{v} \quad \forall \v{v}\times q\in \v{V}_h\times S_h$$ where $\Pi_1$ is the projection from piecewise constant space $P_0$ to continuous piecewise $P_1$, and the constant kernel of the original mixed finite element is gone, however, weird things happened, (2) doesn't work anymore, I coined the test problem from an interface problem for diffusion equation, this is what I got for pressure $p$, the right one is the true solution and the left one is the numerical approximation:

however if $\nu$ is a constant, the test problem performs just fine:

I am guessing it is because the way I am imposing the compatibility condition, since it is linked with the inf-sup stability of the whole system, here is my second question:

• (Q2): is there any way other than (2) to impose the compatibility for pressure $p$? or while coining the test problem, what kind of $p$ should I use?
• MathML not working? – Shuhao Cao Sep 10 '12 at 21:30
• We use MathJaX on StackExchange, everything you posted is showing up beautifully, thanks for the detailed question. – Aron Ahmadia Sep 11 '12 at 8:09

The compatibility condition concerns velocity, not pressure. It states that if you only have Dirichlet boundary conditions for the velocity, then these should be compatible with the divergence-free constraint, i.e. $\int_{\partial \Omega} u \cdot n = 0$ with $\partial \Omega$ the boundary of the computational domain (not the cell).

In this case $\nabla p$ cannot be distinguished from $\nabla (p + c)$ with $c$ an arbitrary constant because you have no boundary condition on $p$ to fix the constant. Thus there are infinitely many solutions for the pressure and in order to compare solutions, a convention is needed. Mathematicians prefer choosing $c$ such that $\overline{p}=p_\mathrm{ref}$ (because they can integrate) while physicist prefer $p(x_\mathrm{ref}) = p_\mathrm{ref}$ (because they can measure in a point). If $Bp$ is your discrete equivalent of $\nabla p$, it implies that $B$ has a null space consisting of the identity vector.

Krylov subspace methods can solve a singular system by removing the null space from the Krylov subspace in which they look for the solution. However, that does not mean you will get the solution $p$ that matches a given convention, you will always need to determine the constant yourself in a postprocessing step, no solver can do it for you.

Here are some suggestions to tackle your problem:

• Equation (2) seems strange. If $\nu$ is a function of $x$ how can it be outside the integral?
• Does your velocity field satisfy the compatibility constraint?
• Try not to do anything for the pressure, just let the solver free to come up with a $p$, then look at $p-p_\mathrm{exact}$. Is it a constant?
• If not, are you sure that the null space of $B$ is indeed the identity vector and nothing more? Both on paper and in the code? The problem seems small enough to actually compute the null space.

As for (Q1), you can choose a solver for saddle-point problems that computes a least squares solution for your system. Then an additional condition can be imposed on the multiplier, like setting a specific degree of freedom, are imposing a specific average.

In general, and I think this answers (Q1), you may use a linear constraint that can distinguish different constants.

This constraint may imposed in a post-processing step, or by a appropiate choice of the trial space (e.g., if you leave out one degree of freedom).