Residual for Stokes equation not equal to zero

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