Problem with solving Stokes equation

Question

Hello,
i try to solve the stoke equation on the unit square and use the the boundary data u=1 on y=1, u=0 on y=0, gN = 0 on x=0 or x=1.
Can you help me to improve my code and to fix the error?

"""-Laplace(u) + grad(p) = f ; div(u) = 0 on the unit square.
u = 1 on y=1.
u = 0 on y=0.
-du/dn - pn = gN 
gN = 0 on x=0 or x=1."""

from dolfin import *

parameters["form_compiler"]["cpp_optimize"] = True
# Create mesh and geometry
mesh = UnitSquareMesh(8,8)

# Define Taylor-Hood function space
V = FunctionSpace(mesh, 'CG', 1)
Q = FunctionSpace(mesh, 'CG', 1)
W= V*Q

# Define function and testfunctions
(u,p) = TrialFunctions(W)
(v,q) = TestFunctions(W)

# Define boundary conditions
tol = 1E-14
def lower_boundary(x, on_boundary):
    return on_boundary and abs(x[1]) < tol

def upper_boundary(x, on_boundary):
    return on_boundary and abs(x[1]-1) < tol

Gamma_0 = DirichletBC(W.sub(0), Constant(0.0), lower_boundary)
Gamma_1 = DirichletBC(W.sub(0), Constant(1.0), upper_boundary)
bcs = [Gamma_0, Gamma_1]


# Define variational problem
f = Constant(0.0)
a = inner(nabla_grad(u), nabla_grad(v))*dx - p*div(v)*dx + div(u)*q*dx 
L = inner(f,v)*dx

# Compute solution
w = Function(W)
solve(a==L,w,bcs)
(u,p) = w.split()

plot(u, titel='Velocity', interactive=True)
plot(p, titel='Pressure', interactive=True)

---

Error message:

BlockquoteAttempting to index with (Ellipsis, Index(13)), but object is already indexed: Indexed(Argument(MixedElement(*[FiniteElement('Lagrange', Domain(Cell('triangle', 2), 'triangle_multiverse', 2, 2), 1, None), FiniteElement('Lagrange', Domain(Cell('triangle', 2), 'triangle_multiverse', 2, 2), 1, None)], **{'value_shape': (2,) }), -2), MultiIndex((FixedIndex(0),), {}))
Traceback (most recent call last):
File "/home/thomas/workspace/Fenics/src/stokes.py", line 50, in
solve(a==L,w,bcs)
File "/usr/lib/python2.7/dist-packages/dolfin/fem/solving.py", line 268, in solve
_solve_varproblem(*args, **kwargs)
File "/usr/lib/python2.7/dist-packages/dolfin/fem/solving.py", line 292, in _solve_varproblem
form_compiler_parameters=form_compiler_parameters)
File "/usr/lib/python2.7/dist-packages/dolfin/fem/solving.py", line 80, in __init__
a = Form(a, form_compiler_parameters=form_compiler_parameters)
File "/usr/lib/python2.7/dist-packages/dolfin/fem/form.py", line 56, in __init__
common_cell)
File "/usr/lib/python2.7/dist-packages/dolfin/compilemodules/jit.py", line 66, in mpi_jit
return local_jit(*args, **kwargs)
File "/usr/lib/python2.7/dist-packages/dolfin/compilemodules/jit.py", line 154, in jit
return jit_compile(form, parameters=p, common_cell=common_cell)
File "/usr/lib/python2.7/dist-packages/ffc/jitcompiler.py", line 77, in jit
return jit_form(ufl_object, parameters, common_cell)
File "/usr/lib/python2.7/dist-packages/ffc/jitcompiler.py", line 159, in jit_form
element_mapping=element_mapping)
File "/usr/lib/python2.7/dist-packages/ufl/form.py", line 123, in compute_form_data
element_mapping=element_mapping)
File "/usr/lib/python2.7/dist-packages/ufl/algorithms/preprocess.py", line 115, in preprocess
form = expand_derivatives(form, common_cell.geometric_dimension()) # EXPR
File "/usr/lib/python2.7/dist-packages/ufl/algorithms/ad.py", line 107, in expand_derivatives
return transform_integrands(form, _expand_derivatives)
File "/usr/lib/python2.7/dist-packages/ufl/algorithms/transformer.py", line 252, in transform_integrands
integrand = transform(integrand)
File "/usr/lib/python2.7/dist-packages/ufl/algorithms/ad.py", line 88, in _expand_derivatives
expression = expand_compounds(expression, dim)
File "/usr/lib/python2.7/dist-packages/ufl/algorithms/expand_compounds.py", line 277, in expand_compounds1
return apply_transformer(e, CompoundExpander(dim))
File "/usr/lib/python2.7/dist-packages/ufl/algorithms/transformer.py", line 277, in apply_transformer
return transform_integrands(e, _transform, domain_type)
File "/usr/lib/python2.7/dist-packages/ufl/algorithms/transformer.py", line 268, in transform_integrands
return transform(expr)
File "/usr/lib/python2.7/dist-packages/ufl/algorithms/transformer.py", line 276, in _transform
return transformer.visit(expr)
File "/usr/lib/python2.7/dist-packages/ufl/algorithms/transformer.py", line 101, in visit
r = h(o, *map(self.visit, o.operands()))
File "/usr/lib/python2.7/dist-packages/ufl/algorithms/transformer.py", line 101, in visit
r = h(o, *map(self.visit, o.operands()))
File "/usr/lib/python2.7/dist-packages/ufl/algorithms/transformer.py", line 101, in visit
r = h(o, *map(self.visit, o.operands()))
File "/usr/lib/python2.7/dist-packages/ufl/algorithms/expand_compounds.py", line 174, in div
return a[...,i].dx(i)
File "/usr/lib/python2.7/dist-packages/ufl/indexed.py", line 94, in __getitem__
error("Attempting to index with %r, but object is already indexed: %r" % (key, self))
File "/usr/lib/python2.7/dist-packages/ufl/log.py", line 154, in error
raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: Attempting to index with (Ellipsis, Index(13)), but object is already indexed: Indexed(Argument(MixedElement(*[FiniteElement('Lagrange', Domain(Cell('triangle', 2), 'triangle_multiverse', 2, 2), 1, None), FiniteElement('Lagrange', Domain(Cell('triangle', 2), 'triangle_multiverse', 2, 2), 1, None)], **{'value_shape': (2,) }), -2), MultiIndex((FixedIndex(0),), {}))

Comments

Please, indent every line of your code by four spaces.

What do you want to do? Stokes problem in domain of topological dimension 2 for scalar velocity and scalar pressure? This is not standard Stokes problem. Does it make sense? Additionally, P1/P1 space is not stable for standard Stokes problem. Is it stable here?

To the errors: UFL is basically complaining that you're taking divergences of a rank 0 (scalar) quantities. Divergence is defined for rank$\ge$1 quantities.

---

At the moment i'm not so familiar with stokes, fem,... I will think about the stabily later.

u should be vector-valued and p scalar?!
=> V = VectorFunctionSpace(mesh, 'CG', 1) ? so i have to fix problem with boundary code

The exact solution of my problem should be given by: u=(y,0); p=0

Answer 1

Jan is right - velocity must come from a vector space. Change

V = FunctionSpace(mesh, 'CG', 1)

to

V = VectorFunctionSpace(mesh, 'CG', 1)

Also, look at the DOLFIN Stokes demos.

Comments

I adapt the part of the boundary condition and now it works.
But i'm not sure if i did it correct?! (u=1 on y=1, u=0 on y=0, gN = 0 on x=0 or x=1).

def lower_boundary(x, on_boundary):
    return on_boundary and abs(x[1]) < tol

def upper_boundary(x, on_boundary):
    return on_boundary and abs(x[1]-1) < tol

Gamma_0 = DirichletBC(W.sub(0), (0.0,0.0), lower_boundary)
Gamma_1 = DirichletBC(W.sub(0), (1.0,0.0), upper_boundary)
bcs = [Gamma_0, Gamma_1]

---

As result i get a function u. And I want to compare the solution with the exact result u=(y,0). How I can do that?