Residual for Stokes equation not equal to zero

1

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
2

I forgot to add that I use docker and my FEniCS version is 2019.1.0.

3

I’m not convinced you’re computing the residual. Cf.

# Print residual
r = assemble(a(u, p, v, q) - L(v, q))
for bc in boundaries:
    bc.apply(r, w.vector())
print(r.norm("l2"))

>>> 9.950073998023825e-14
2 Likes
4

Thank you very much for the answer! I thought as well that it’s something about the boundary conditions.

However, I have a follow-up question. In this example, as far I understood, this is a residual over the whole test space. What I’m interested in, is computing the residual for a specific function. What should I do to correctly compute the residual for an arbitrary function? I took (u, p) as an example since from the definition of the problem we have that (u, p) \in W and thus I should get a((u, p)(u, p)) - L(u, p) = 0. Using

for bc in boundaries:
    bc.apply(r, w.vector())

doesn’t work, because then I need to have r defined using GenericVector or a GenericMatrix and not a “normal” Function. Are there any other ways to apply boundary conditions?

5

I’m a bit confused regarding exactly what you want to achieve. Could you clarify what you want to measure? What you’ve written looks close to the energy norm of the true solution and not the residual.

6

So, I think @Martyna_Soszynska’s mistake was to assume that w is in the test space. The inflow BC is nonzero, but, even when Dirichlet data is nonzero, the test functions still go to zero at the Dirichlet boundaries. Thus, the solution satisfying the inflow BC is outside the space of test functions for which the residual is zero. If we instead apply homogenized versions of the BCs to a copy of w, then plug that in as a test function, we get the expected result:

# ...

# Compute solution
w = Function(W)
solve(a(u, p, v, q) == L(v, q), w, boundaries)

# Split the mixed solution
(u, p) = w.split()

# Get a test function with homogeneous versions
# of all BCs; even when the solution satisfies a
# nonzero Dirichlet BC, the test space goes to
# zero at the Dirichlet boundary.
w0 = Function(W)
w0.assign(w)
for b in boundaries:
    b.homogenize() # (Sets Dirichlet data to 0)
    b.apply(w0.vector())
(v,q) = w0.split()
    
# Print residual
#print(assemble(a(u, p, u, p) - L(u, p)))
print(assemble(a(u, p, v, q) - L(v, q)))
1 Like
7

@kamensky, thank you, I think this is exactly what I was missing :slight_smile:

So the thing that I didn’t get is that the trial and the test space are not the same - even if I set nonzero Dirichlet boundary conditions on the trial space, the test space will still have zero on the boundary. So whenever I want to take something from the trial space and put it in the position of the test space (or vice versa), I need to take it into account. In this case, when I want to take the solution and put it in the position of the test function, I need to set the boundary conditions to zero.

@nate, I work with so-called Dual Weighted Residual Method for an a posteriori error estimation. There, I take a functional J(u, p), where (u, p) is the exact solution of the problem and want to estimate |J(u, p) - J(u_h, p_h)|. The idea is to take

\rho((u_h, p_h)(v_h, q_h)):= a((u_h, p_h)(v_h, q_h)) - L(v_h, q_h)

and instead of (v,_h q_h) as the second argument plug in some other well chosen function. Actually it is (z^{(1)}_h - z_h, r^{(1)}_h - r_h), where (z_h, r_h) is the solution of the discrete adjoint problem, and (z^{(1)}_h, r^{(1)}_h) is the reconstruction of the exact adjoint solution. Then it will estimate the error |J(u, p) - J(u_h, p_h)|. In this method, \rho is called primal residual (the adjoint counterpart is called adjoint reasidual). So that’s why I used the term residual in this post.

I didn’t want to put my lengthy code, so that’s why I chose the Stokes example to illustrate my problem. But as a trade-off there computing a((u_h, p_h)(u_h, p_h)) - L(u_h, p_h) didn’t make much sense and it could have been confusing.

It’s been puzzling me for long. In my program I need to juggle between functions coming from trial and test spaces. And as soon as I took a problem with nonzero Dirichlet boundary conditions the results were off and I had no idea why. So thank you very much, both of you :slight_smile:

8

I have implemented the dual problem computation and goal oriented DWR estimators for the Poisson, Burgers and Compressible Navier Stokes systems. Perhaps they’ll help. Note that I enrich the dual space by increasing the polynomial order, which is the expensive and lazy option.