Residual for Stokes equation not equal to zero

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