My basic Stokes Example won't work.

Question

This code produces a solution that doesn't match the example I cooked up.

My exact solution is: u = (sin(pix)cos(piy),sin(piy)cos(pix))
My exact pressure is: p = picos(piy)cos(pix)

I think I'm imposing a boundary condition incorrectly, because when I plot "my" solutions, the boundary values are way off.

Any tips would be of great help. Thank you so much!

from dolfin import *
mesh = UnitSquareMesh(32,32)
#Function Spaces
V = VectorFunctionSpace(mesh, "Lagrange", 2)
Q = FunctionSpace(mesh, "Lagrange", 1)
W = V*Q


#Trial and Test Functions
(u,p) = TrialFunctions(W)
(v,q) = TestFunctions(W)

#Exact Values / RHS values
f = Expression(("-pi*pi*sin(pi*x[0])*cos(pi*x[1])","-pi*pi*sin(pi*x[1])*cos(pi*x[0])"))
u_ex = Expression(("sin(pi*x[0])*cos(pi*x[1])","sin(pi*x[1])*cos(pi*x[0])"))
p_ex = Expression("pi*cos(pi*x[0])*cos(pi*x[1])")

#Boundary Conditions
def right(x, on_boundary): return x[0] > 1.0 - DOLFIN_EPS
def left(x, on_boundary): return x[0] < DOLFIN_EPS
def top(x, on_boundary): return x[1] > 1.0 - DOLFIN_EPS
def bottom(x, on_boundary): return x[1] < DOLFIN_EPS    

bc_u_right = DirichletBC(W.sub(0), u_ex, right)
bc_u_left = DirichletBC(W.sub(0), u_ex, left)
bc_u_top = DirichletBC(W.sub(0), u_ex, top)
bc_u_bottom = DirichletBC(W.sub(0), u_ex, bottom)

bc_p_right = DirichletBC(W.sub(1), p_ex, right)
bc_p_left = DirichletBC(W.sub(1), p_ex, left)
bc_p_top = DirichletBC(W.sub(1), p_ex, top)
bc_p_bottom = DirichletBC(W.sub(1), p_ex, bottom)

bcs = [bc_u_right, bc_u_left, bc_u_top, bc_u_bottom, bc_p_right, bc_p_left, bc_p_top, bc_p_bottom]


#Variational Problem
a = inner(nabla_grad(u), nabla_grad(v))*dx - div(v)*p*dx + q*div(u)*dx
L = inner(f, v)*dx

#Solve
w = Function(W)
solve(a == L,w)
(u,p) = w.split()

#Plot
plot(u, title="My Velocity")
plot(u_ex,mesh, title="Exact Velocity")

plot(p, title="My Pressure")
plot(p_ex, mesh, title="Exact Pressure")

interactive()

Answer 1

Hi, first of all your exact velocity is not divergence free, so you are not testing against an analytical solution of the Stokes problem. Further, you must not set both velocity and pressure on the same part of boundary. Finally, you forgot to pass the list of boundary conditions to the solve function. Here's a modified example

from dolfin import *
mesh = UnitSquareMesh(32,32)
#Function Spaces
V = VectorFunctionSpace(mesh, "Lagrange", 2)
Q = FunctionSpace(mesh, "Lagrange", 1)
W = V*Q

#Trial and Test Functions
(u,p) = TrialFunctions(W)
(v,q) = TestFunctions(W)

#Exact Values / RHS values

u_ex = Expression(('cos(pi*x[1])', 'sin(pi*x[0])'))
p_ex = Expression('pi*cos(pi*x[0])*cos(pi*x[1])')

f = Expression(('-pi*pi*sin(pi*x[0])*cos(pi*x[1]) + pi*pi*cos(pi*x[1])',
                'pi*pi*sin(pi*x[0]) - pi*pi*sin(pi*x[1])*cos(pi*x[0])'))

#Boundary Conditions
def right(x, on_boundary): return x[0] > 1.0 - DOLFIN_EPS
def left(x, on_boundary): return x[0] < DOLFIN_EPS
def top(x, on_boundary): return x[1] > 1.0 - DOLFIN_EPS
def bottom(x, on_boundary): return x[1] < DOLFIN_EPS

bc_u_right = DirichletBC(W.sub(0), u_ex, right)
bc_u_left = DirichletBC(W.sub(0), u_ex, left)
bc_u_top = DirichletBC(W.sub(0), u_ex, top)
bc_u_bot = DirichletBC(W.sub(0), u_ex, bottom)

bcs = [bc_u_right, bc_u_left, bc_u_top, bc_u_bot]

#Variational Problem
a = inner(nabla_grad(u), nabla_grad(v))*dx - div(v)*p*dx + q*div(u)*dx
L = inner(f, v)*dx

#Solve
w = Function(W)
solve(a == L, w, bcs)
(u, p) = w.split(True)
normalize(p.vector())

#Plot
plot(u, title="My Velocity")
plot(u_ex,mesh, title="Exact Velocity")

plot(p, title="My Pressure")
plot(p_ex, mesh, title="Exact Pressure")

interactive() 

Note that since all bcs are set on velocity, pressure is only determined up to a constant. The solution is made unique by normalization. You might want to check the Stokes demo.