Optimizations in 'quadrature' mode are not safe w.r.t. floating point rounding errors

[Martin Sandve Alnæs]

from dolfin import *

if 0: # Works
    parameters["form_compiler"]["representation"] = "uflacs"
if 0: # Works
    parameters["form_compiler"]["optimize"] = False
    parameters["form_compiler"]["representation"] = "quadrature"
if 1: # Fails
    parameters["form_compiler"]["optimize"] = True
    parameters["form_compiler"]["representation"] = "quadrature"

n = 128
mesh = UnitSquareMesh(n, n)

class Source(Expression):
    def eval(self, values, x):
        values[0] = 4.0*pi**4*sin(pi*x[0])*sin(pi*x[1])

x = SpatialCoordinate(mesh)
uexact = sin(pi*x[0])*sin(pi*x[1])
f1 = div(grad(div(grad(uexact))))
f2 = 4.0*pi**4*sin(pi*x[0])*sin(pi*x[1])
f3 = Source(degree=2)

e13 = assemble((f1 - f3)**2*dx)
e23 = assemble((f2 - f3)**2*dx)
e12 = assemble((f1 - f2)**2*dx)

print("Expecting the three following numbers to be ~0:")
print(e13)
print(e23)
print(e12)
assert e13 < 1e-8
assert e23 < 1e-8
assert e12 < 1e-8  # Fails for optimized quadrature, ~4000 instead of ~0