Hazardous floating point simplifications

[Jan Blechta]

"Small" floats are dropped in FFC so that the following code

coord = VectorElement("P", triangle, 1)
mesh = Mesh(coord)
dx = dx(mesh)

M = 1e-16*dx

results in

    // Compute element tensor
    A[0] = 0.0;

1e-16 is perfectly good floating point number. It is negligible compared to 1.0 but the compiler should not assume that it will be compared to 1.0.

Recommended reading: What every computer scientist should know about floating-point arithmetic (Goldberg, 1991).

[Martin Sandve Alnæs]

I guess it still makes sense to clamp values to zero in element tables, where we know the source of the values, but not for arbitrary floats.

Try adjusting the precision/epsilon parameters. Maybe there should be one precision for tables and one for other floats.

Note that one of the reasons for lower precision is reproducable sources and results in the regression tests, this must be preserved.

[Jan Blechta]

I agree. Numerical zeros in element tables which should be exact zeros should better be removed. I guess that this numerical artifact comes from FIAT loosing some structure due to traditional reference element coordinates.

I believe that setting smaller epsilon will help. But that's a workaround which is pretty much hidden from user expecting to integrate the form above exactly. We need to check the smallness tests in FFC for correctness.

I know about the testing issue. In fact, the lack of precision parameter in tsfc blocks us passing the tests. @miklos1, the default values of EPSILON in constants.py and epsilon in fem.py are good. But could those be overridden by parameters precision and epsilon? (Note that PRECISION is used in coffee.py and epsilon in fem.py.) It seems that monkey-patching PRECISION did exactly what we need.

[Miklós Homolya]

I haven't checked, but I believe the example in the issue description works with TSFC, since rounding to zero is done for element tables only. Unfortunately, this isn't so on the FInAT branch at the moment, since FInAT returns expressions instead of tabulation matrices, so it becomes less obvious which numbers should be rounded to zero, so all are rounded for now.

I think the ideal solution would be to make FIAT precise (say, up to machine precision). sympy has arbitrary precision floating-point numbers, although I'm not sure if one can only prescribe the precision of individual floating-point operations, or also the precision of the result of the numerical evaluation of symbolic expressions. Evaluating the entries of element tables up to the precision of double should yield exactly zero in all cases when the symbolic expression could simplify to zero (even if sympy cannot figure out the simplification symbolically).

default values of EPSILON in constants.py and epsilon in fem.py are good.

Yes, they were chosen to match FFC's default values because without these floating-point simplifications some Firedrake tests were failing due slightly higher numerical error than the prescribed numerical tolerance.

But could those be overridden by parameters precision and epsilon?

Can you file an issue for this on GitHub?

[Jan Blechta]

@miklos1, https://github.com/firedrakeproject/tsfc/issues/61.

[Jan Blechta]

NOTE for me: Miklos summarized rounding procedures in TSFC/COFFEE https://github.com/firedrakeproject/tsfc/pull/66#issuecomment-253770891. See also https://github.com/firedrakeproject/tsfc/pull/76.

[Martin Sandve Alnæs]

I've taken a look at uses of tolerances in uflacs in a branch soon to be merged, here'a a decent manual test:

In [1]: from dolfin import *
   ...: #parameters["form_compiler"]["representation"] = "uflacs"
   ...: mesh = UnitSquareMesh(16, 16)
   ...: c = Constant(1e+16)
   ...: M = (1.0 - c*1e-16)*dx(mesh)
   ...: print(assemble(M))
   ...: 
Calling FFC just-in-time (JIT) compiler, this may take some time.
Calling FFC just-in-time (JIT) compiler, this may take some time.
Calling FFC just-in-time (JIT) compiler, this may take some time.
Calling FFC just-in-time (JIT) compiler, this may take some time.
1.0

vs

In [2]: from dolfin import *
   ...: parameters["form_compiler"]["representation"] = "uflacs"
   ...: mesh = UnitSquareMesh(16, 16)
   ...: c = Constant(1e+16)
   ...: M = (1.0 - c*1e-16)*dx(mesh)
   ...: print(assemble(M))
   ...: 
Calling FFC just-in-time (JIT) compiler, this may take some time.
Calling FFC just-in-time (JIT) compiler, this may take some time.
Calling FFC just-in-time (JIT) compiler, this may take some time.
Calling FFC just-in-time (JIT) compiler, this may take some time.
0.0

[Jan Blechta]

Possibly WONTFIX as uflacs will be adopted as a default.

[Jan Blechta]

UFLACS and TSFC do not fail on the MWE. Quadrature does, but it will be hopefully removed soon.