Hello everyone,

I work on error estimation based on residual methods. Thus, unsurprisingly, it’s important for me to compute residuals accurately. And I have a problem with Stokes equation.

In a variational problem we have something like:

Find (u, p) \in W such that:

a((u, p)(v, q)) = L(v, q)

for all (v, q) \in W

So, then what I should get after solving the problem is a((u, p)(u, p)) - L(u, p) = 0 (or something very close to 0). But in the example below I get something much larger, namely 1.3333333333337636. I Why is it? I have a feeling that I’m missing something.

Here is an easy sample of my problem (based on the tutorial https://fenicsproject.org/docs/dolfin/1.4.0/python/demo/documented/stokes-taylor-hood/python/documentation.html)

```
from dolfin import *
# Load mesh and subdomains
mesh = RectangleMesh(
Point(0.0, 0.0), Point(4.0, 1.0), 80, 20, diagonal="right"
)
# Define boundaries
def boundary_down(x, on_boundary):
return on_boundary and near(x[1], 0.0)
def boundary_up(x, on_boundary):
return on_boundary and near(x[1], 1.0)
def boundary_left(x, on_boundary):
return on_boundary and near(x[0], 0.0)
def boundary_right(x, on_boundary):
return on_boundary and near(x[0], 4.0)
# Define function spaces
V = VectorElement("Lagrange", mesh.ufl_cell(), 2)
Q = FiniteElement("Lagrange", mesh.ufl_cell(), 1)
W = FunctionSpace(mesh, V * Q)
# Define Dirichlet boundaries
inflow = Expression(("x[1] * (1.0 - x[1])", "0.0"), degree=2)
boundaries = [
DirichletBC(W.sub(0), Constant((0.0, 0.0)), boundary_up),
DirichletBC(W.sub(0), inflow, boundary_left),
DirichletBC(W.sub(0), Constant((0.0, 0.0)), boundary_down),
DirichletBC(W.sub(1), Constant(0.0), boundary_right),
]
# Define variational problem
(u, p) = TrialFunctions(W)
(v, q) = TestFunctions(W)
f = Constant((0, 0))
# Define forms
def a(u, p, v, q):
return (inner(grad(u), grad(v)) - div(v) * p + q * div(u)) * dx
def L(v, q):
return inner(f, v) * dx
# Compute solution
w = Function(W)
solve(a(u, p, v, q) == L(v, q), w, boundaries)
# Split the mixed solution
(u, p) = w.split()
# Print residual
print(assemble(a(u, p, u, p) - L(u, p)))
# Save solution in VTK format
ufile_pvd = File("velocity.pvd")
ufile_pvd << u
pfile_pvd = File("pressure.pvd")
pfile_pvd << p
```