Computational Science Stack Exchange is a question and answer site for scientists using computers to solve scientific problems. It's 100% free, no registration required.

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

I am looking for an easy to implement stabilization scheme that can be used with equal order ($P_1-P_1$ or $Q_1-Q_1$) finite elements for fluid flow. Is there something like this or should I stick to Taylor Hood elements? It appears that there are a number of methods but the ease of implementation is not quite clear.

share|improve this question


First, be warned that stabilization is generally much more fragile than use of stable elements. If you solve a narrow class of problems, you might be able to choose reasonable parameters and get on with your work, but if you solve a wide range of problems or need to solve demanding under-resolved scenarios, especially with less regular boundaries, coefficient structure, and constitutive relations, stabilization will generally require constant vigilance to detect and quantify the impact of numerical artifacts, and to re-tune parameters for.

For these reasons, I strongly recommend stable elements (especially those with discontinuous pressure space, thus locally conservative) if you want a robust method that does not need to be frequently revisited. In my opinion, reliability and clean analysis is worth the modest increase in data structure complexity.

Stabilization methods

I'll outline the two most popular choices in the context of a stabilization term $S_h(\mathbf u_h, p_h, \mathbf f; \mathbf v_h, q_h)$ to be added to the weak form.

  • Pressure-stabilized Petrov-Galerkin (PSPG) This stabilization has the form $$ S_h^{\text{pspg}}(\mathbf u_h,p_h;q_h) = \sum_{T \in \mathcal{T}_h} \alpha_T\Big(-\nabla\cdot(\eta D \mathbf u_h) + \nabla p_h - \mathbf f, \nabla q_h\Big)_T$$ where the inner product $(\cdot,\cdot)_T$ is evaluated over element $T$. This method was introduced by Hughes, Franca, and Balestra in 1986 and is a residual-based method. It requires second derivatives on elements (if applicable; they are only zero for lowest order spaces with affine elements). When the solution is not very well resolved, the element integral can be very different from the actual cell residual (because the influence from neighbor cells is discarded). The stabilization produces a non-physical boundary condition for the pressure, leading to local inaccuracies impacting functionals of interest (such as drag). There are more robust and accurate edge-stabilized modifications, but these are more complicated and grow the stencil. Residual stabilization adds a large number of couplings for more complicated problems (like reaction and turbulence models) and significant implementation complexity. The formulation is cumbersome for time-dependent problems because the "affine" term enters in an inner product, thus typically precluding a method of lines approach. PSPG also breaks symmetry (relevant for adjoint-consistency) and increases regularity requirements for an adjoint, making it less attractive for optimization and sensitivity analysis.

  • Local pressure projection This conceptually simpler method uses a stabilization of the form $$ S_h^{\text{lps}} = (\alpha \tilde \pi p_h, \tilde \pi q_h)$$ where $\tilde \pi$ is a local fluctuation operator and $\alpha$ is a scaling factor. Dohrmann and Bochev (2004) studied what is perhaps the simplest choice, $\tilde \pi = 1 - \pi_{k-1}$, where $\pi_{k-1}$ is cell-wise projection into a polynomial space of degree one less than the continuous space. There are many patch-based and edge-stabilized choices, see Braack and Lube's 2009 review. These methods are simpler to implement and preserve symmetry, but are often less robust.

share|improve this answer
Thanks for the detailed answer. Yes the class of problems is rather narrow. I will look more into the pressure projection method. At some point the code will be used for inversions so being able to solve the adjoint problem will be useful. – stali Jan 17 '13 at 23:17
I looked at some papers but how exactly the stabilization matrix is constructed is not clear to me. – stali Feb 17 '13 at 23:43
@stali There are many papers, books, and open source implementations available. If after exhausting these resources, you are still stuck on something, you're welcome to ask more focused questions on this site. – Jed Brown Feb 18 '13 at 0:19

This is a late response, but I want to point out the option (for the time-dependent case) of using a pressure-approximation type of projection method, based on solving a Poisson equation for the pressure with a well-posed boundary condition. See [1] for derivations and tests, including tests with equal-order finite elements.

I have to admit that Jed Brown's warning still largely applies, however --- the schemes proposed and discussed in [1] appear very flexible and are supported by some theory and some testing, but have not been widely tested on many different problems. They may prove to be robust, but this has to be proved by experience.

[1] J.-G. Liu, J. Liu, R.L. Pego, Stable and accurate pressure approximation for unsteady incompressible viscous flow, J. Comp. Phys. 229 (2010) 3428-3453.

share|improve this answer

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

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