Solving Stokes equation with piecewise constant pressure

Question

Looking at the online Stokes examples, (e.g., number 16 with iterative solver), I want to apply P2/P0 and perhaps P2/(P1+P0) element combination. For instance, in the latter case I figured all I had to do was change this

V = VectorFunctionSpace(mesh, "CG", 2)    
Q = FunctionSpace(mesh, "CG", 1)
W = V * Q

to

V = VectorFunctionSpace(mesh, "CG", 2)  
Q = FunctionSpace(mesh, "CG", 1)  
P = FunctionSpace(mesh, "DG", 0)
W = V * (Q + P)

While leaving everything else the same. For P2/P0 I would simply change Q to "DG" and order 0. But when I do, I get the following pressure results (LHS is P2/(P1+P0), RHS is P2/P1 )

If I did P2/P0 I also get a similar result as the LHS image. My velocity solutions are not negatively affected (in fact, they are "better" because they are now locally conservative).
How would I solve a P2/(P1+P0) (or even a P2/P0) such that my pressure solutions do not get screwed up? Would it have something to do with how I apply my pressure boundary conditions?

Comments

Can you provide more detail why you need local conservation properties? Do you want to transport something with your flow field?

Answer 1

The P2-P0 is not LBB stable in 3D, and thus also the enriched Taylor-Hood cannot be stable. This may explain why your pressure does not look as expected. In fact, achieving local conservation in 3D is not as straightforward as in 2D, since most tetrahedral element pairings with discontinuous pressures are notoriously unstable. You may try P2+Bubble paired with DG1 though (see Crouzeix-Raviart 1973 or so). However I have never encountered this element in the wild since it is pretty expensive to deal with in most general purpose codes. If you need local conservation you may also try Hdiv-spaces with DG-stabilization for the velocity (e.g. as in some of Guido Kanschat's works). This works straightforwardly in FEniCS I would guess. However this also may be tough for blackbox linear solvers.