Difference in computed residuals

1

Hi everyone, I need to check some residuals but I didn’t see convergence. After some tests, I figured that there is a strange behavior, that I report below. Perhaps it is stupid, but thanks for the help.

from dolfin import *

Mesh

N = 10
mesh = UnitSquareMesh(N, N)

Functions

V = FunctionSpace(mesh, ‘CG’, 1)
bc = DirichletBC(V, Constant(0), ‘on_boundary’)

u = Function(V)
v = TestFunction(V)

f = Constant(1)

Residual and matrices

aa = dot(grad(u), grad(v))
ff = f*v
F = (aa - ff)*dx
dF = derivative(F, u, TrialFunction(V))
A_res = assemble(dF)
b_res = assemble(-F)

aa = dot(grad(TrialFunction(V)), grad(v))
A = assemble(aadx)
b = assemble(ffdx)
for obj in (A, b, A_res, b_res):
bc.apply(obj)

sol_res = Function(V)
sol = Function(V)

solve(A_res, sol_res.vector(), b_res)
solve(A, sol.vector(), b)

Compute residual norms

Fvec = PETScVector()
assemble(F, tensor=Fvec)
print(“Difference:”, (sol.vector() - sol_res.vector()).vec().norm()) # 0
print(“Residual F”, Fvec.vec().norm()) # 0.09
print(“Residual Ax-b”, (A*sol.vector()-b).vec().norm()) # 1e-16

Thanks in advance

2

When you assemble the Form F at the end, it is using the Function u which was used to define it. F is not connected in any way to either sol or sol_res. If you run

u.assign(sol_res)

immediately before assembling F, then apply the BC (or, more generally, a homogenized version of the BC, if it were not zero) to Fvec, you get that the norm of Fvec is machine zero.

3

thanks! That was careless. Indeed, it works this way. Still, norms aren’t exactly equal, but now it works.

Thanks again,
Nico