Laplacian in the stabilizatio term of Badia & Codina

Question

I'm trying to implement the stabilization techniques of Badia & Codina, based on the variational multiscale concept. For some equations, in particular the Stokes problem, this
technique leads to stabilization terms such as:

laplace = lambda v: grad(div(v)) - curl(curl(v))

a = .... + (-mu*laplace(u) + grad(p), mu*laplace(v) + grad(q))*dx + ...

where a is the bilinear form, u,p are trial functions for velocity and pressure, v, q are the
corresponding test functions.

The method fails spetacularly whenever the laplacian in the bilinear form is actually invoked, ie, with elements of order higher than 1.

Amazingly enough, the stabilization seems to work even if I remove these laplacian terms, or even (more amazingly still) if I remove only the laplacian term of the trial function or even (more more amazingly still) if I replace the laplacian operator on the velocity with an indentity function. I'm not proposing these are general results :)

My question is: is there a known issue with using the second order derivatives in the bilinear form? In a similar question, someone suggested using integration by parts, but this won't help me because there is a product of laplacians, I could perhaps just disregard this product as a higher order term, and do the integration by parts, but I'd like to hear some opinions on the whole situation.

Comments

Dear rbw,

please could you provide a minimal working example producing the failure (full code of it)?

Regards,
Peter.

---

The following example shows the situation.
Just for reference, I'm trying to implement:
Badia, Codina, Unified Stabilized Finite Element Formulations for the Stokes and the Darcy Problems, SIAM, 2008.
This is a Stokes problem, but for a Darcy problem the same behaviour happens.
You just have to replace the statement that sets the order of the elements and look at the solution.

 # Half side length
len=0.50
mesh = RectangleMesh(-len, -len, len, len, 30, 30) 
mu = 1.0

order = 1 # Ok, even if not super convergent :)
# order = 3 # Catastrophe, but should work
V = VectorFunctionSpace(mesh, "CG", order)
Q = FunctionSpace(mesh, "CG", order)
W = V * Q

# Exact Solutions
velocity = Expression(("pi*sin(4*pi*x[0])*cos(4*pi*x[1])",
                      "-pi*cos(4*pi*x[0])*sin(4*pi*x[1])"))
pressure = Expression("8*pi*pi*sin(4*pi*x[0])*cos(4*pi*x[1])")
# f = -mu lambda u + grad p
f = Expression(("32*(pi*pi*pi)*(cos(4*pi*x[1])*sin(4*pi*x[0]) + \
                                cos(4*pi*x[0])*sin(4*pi*x[1]))", 
                "32*(pi*pi*pi)*(cos(4*pi*x[1])*sin(4*pi*x[0]) - \
                                cos(4*pi*x[0])*sin(4*pi*x[1]))"))

class NoSlipDomain(SubDomain):
  def inside(self, x, on_boundary):
    return on_boundary

bc0 = DirichletBC(W.sub(0), velocity, NoSlipDomain())

# Trial and test functions
(u, p) = TrialFunctions(W)
(v, q) = TestFunctions(W)
w = Function(W)

# Stabilization parameters
h = CellSize(mesh)
c1 = 1.0
taup = c1*mu
tauu = (h*h)/(c1*mu)

# Define the variational problem
# laplace = lambda v: div(grad(v))
laplace = lambda v: grad(div(v)) - curl(curl(v))

a = mu*inner(grad(u),grad(v))*dx - p*div(v)*dx + q*div(u)*dx \
  + taup*div(v)*div(u)*dx \
  + tauu*inner(grad(q), grad(p))*dx \
  + tauu*inner(mu*laplace(v)+grad(q), -mu*laplace(u)+grad(p))*dx
# + jump terms

L = inner(f, v)*dx \
  + tauu*inner(f, mu*laplace(v)+grad(q))*dx
# + jump terms

A, b = assemble_system(a, L, bc0)
solve(A, w.vector(), b)
u, p = w.split(True)
p.vector()[:] -= p.vector().array().mean()

plot(u)
plot(p)
plot(interpolate(pressure, Q))
interactive()

---

Hello,

I tried to solve this problem and realized that for pressure in space of order order and velocity in space of order order - 1 the results are "quite" alright.

But you are right, there must be probably a glitch in a numerical integration. But I am just a newbie and here are more erudite persons, hope somebody else will tell you more.

I also tried higher orders (again with plus one higher order for pressure) and it worked, as far as I just see the results.

---

I upvoted this question at least.

Answer 1

Perhaps your stabilisation parameters should depend on the polynomial order? The results looks much better with cubic elements if tauu is is multiplied by 3.