nan values when computing arccos(1.0) - bug?

Question

This code is pretty self-explanatory:

from dolfin import *
import numpy as np
mesh = UnitSquareMesh(1, 1)
S = FunctionSpace(mesh, 'CG', 1)
r = interpolate(Expression('1.0', degree=1), S)

# with fenics: result is [ nan  nan  nan  nan]
acos_r = acos(r)
print project(acos_r,S).vector().array()

# with numpy: result is [ 0.  0.  0.  0.]
print np.arccos(r.vector().array())

I assume there is some numerical issue since the input to arccos has an upper bound of <= 1 (in contrast, using e.g. 0.999 returns the same answer for both methods). Is there a simple workaround?

Answer 1

Indeed, probably a numerical issue. But can't you just get away with adding a small 'eps' as:

r = interpolate(Expression('1.0-DOLFIN_EPS', degree=1), S)

Which leads to a result numerically zero.

Now, if you want to bound it both from above and from below, you can consider the following approach:

from dolfin import *
import numpy as np
mesh = UnitSquareMesh(1,1)
S = FunctionSpace(mesh, 'CG', 1)

# Construct (conditional) expression
DOLFIN_EPS = 1E-14
v    = Expression('cos(2*pi*x[0])',degree = 1)
# Augment with conditionals to restrict solution from above and from below:
c1  = conditional( ge(v,  1-DOLFIN_EPS), 1-DOLFIN_EPS, v )
c2  = conditional( le(v, -1+DOLFIN_EPS), -1+DOLFIN_EPS, c1 )

r = project(c2, S)

# Now you get finite results, both in FEniCS and with numpy!
acos_r = acos(r)
print project(acos_r,S).vector().array()
print np.arccos(r.vector().array())

Note that you actually can skip one projection step, so in the above presented code you can also do it just like this:

 acos_r = acos(c2)
 print project(acos_r,S).vector().array()

Furthermore, note that it really is a workaround. So no guarantees on robustness ;)

Comments

For an arbitrary function in the range [-1, 1], do you know if there's a simple way to threshold the function e.g. between [-1+DOLFIN_EPS, 1-DOLFIN_EPS]?

---

See the extended answer. Hope this helps!

---

Thanks! Haven't used conditionals before so this was really useful for me :)