Checking derivative of a form in the presence of boundary conditions

Question

I'm computing the Jacobian of a form using derivative() and I'd like to double-check the result using simple finite differences. Here's my setup:

mesh = UnitSquareMesh(2, 2)
U = FunctionSpace(mesh, "Lagrange", degree=1)
V = U * U
uf = Function(V)
u, f = split(uf)
v, w = split(TestFunction(V))
a = dot(grad(v), grad(u)) * dx - v * f * dx
da = derivative(a, uf)

The values obtained by assembling da are consistent in the sense that they're close to values obtained by applying simple finite differences to a.

Now if I introduce boundary conditions

u0 = interpolate(Constant(0), U)
bc = DirichletBC(U, u0, lambda x, on_boundary: on_boundary)

and check da at a point where I apply the boundary conditions, the values aren't consistent any more. For example:

import numpy as np
x = np.ones_like(uf.vector().array())
uf.vector()[:] = x
bc.apply(uf.vector())  # Apply bc to uf
cx = assemble(a).array()
Jx = assemble(da).array()

e = np.zeros_like(x); e[0] = 1
h = 1.0e-6

uf.vector()[:] = x + h * e
bc.apply(uf.vector())
cxph = assemble(a).array()

uf.vector()[:] = x - h * e
bc.apply(uf.vector())
cxmh = assemble(a).array()

diff = (cxph - cxmh) / (2 * h)
print diff[0]
print Jx[0, 0]

This code outputs:

0.0
1.0

How can I ensure that I only evaluate a and da at values satisfying the boundary conditions and have consistent derivatives?

Thanks!

Answer 1

The finite differences gives you an approximate directional derivative (a vector), while Jx is the derivative operator (a matrix). In order to compare the two you need to apply the matrix to the same directional vector you used to compute the finite differences. In your example this vector is a zero-vector.

Also note that the correct way to handle Dirichlet boundary conditions on subspaces of a mixed function space is

bc = DirichletBC(V.sub(i), u0, subdomain)

where i is the index of the subpace.

Comments

Thanks. In my code, diff[0] is an approximation to the derivative of the "first" constraint with respect to the "first" variable, which is what Jx[0, 0] should be, I believe?! If I don't apply the boundary conditions, they match very closely.

Thanks for the tip on BCs in mixed space.