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

Question

$\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?

Answer 1

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.

Answer 2

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