floating-point based issues for replays in firedrake

Hi all,

tests_firedrake/poisson_bc/ currently fails (sometimes?!) on mac systems.

The reason appears to be the following.

Firedrake always solves variational problems in residual form (lifting the boundary conditions to the right hand side and solving a homogeneous problem).

To match up with this for solves with assembled operators, we basically do the same. But because we've been given an assembled right hand side, we have to lift the boundary values at the linear algebra level.

This leads to floating point mismatches when using replay_dolfin with a tol=0 (as used in the referenced test). It's demonstrated by this example:

from firedrake import *
mesh = UnitSquareMesh(1, 1)
V = FunctionSpace(mesh, "CG", 1)

u = TrialFunction(V)
v = TestFunction(V)
a = u*v*dx
L = v * dx

u = Function(V, name="u")

# Linear operators
f = assemble(assemble(L) - assemble(action(a, u)))

# Variational solve
g = assemble(L - action(a, u))

print f.dat.data - g.dat.data

Effectively, libadjoint (implicitly) assumes that floating point is associative in the face of such manipulations. This is sort of unsafe.

I wonder if the thing to do is just to relax the tolerance in these tests.

I note in passing that the equivalent code in dolfin, when using the uflacs representation, also suffers from this "issue".

Comments (2)