Divide by zero error when using conditionals

[Tanyakarn Treeratanaphitak]

Hi,

I'm trying to solve a nonlinear problem that requires me to compute the value of 1/u, where u is my unknown. In the event that u is near zero, I used a conditional statement to set the value of that particular term to zero. However, when I tried running my simulation, I'm getting the residuals to be NaNs.

Here is an example of what I did:

norm = sqrt(dot(u,u))
f = conditional(le(norm, 0.001), 0., conditional(And(lt(norm, 10.), gt(norm, 0.001)), (2./norm)*(1. +  pow(norm,0.7)), 0.14)) * u

F = inner((1 + u**2)*grad(u), grad(v))*dx - f*v*dx

Output from FEniCS:

Solving nonlinear variational problem.
Newton iteration 0: r (abs) = 5.745e+00 (tol = 1.000e-10) r (rel) = 1.000e+00 (tol = 1.000e-06)
Newton iteration 1: r (abs) = -nan (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-06)

I took a look at the generated header file and it seems that what FEniCS is doing is (-(conditional result as 0/1) + 1) * (expression) and since I have something divided by zero, it does not matter what the first part returns, the result will still be NaN. I get that it's to avoid evaluating the conditional, but when the conditional there to avoid another issue, it needs to be evaluated.

Thank you.

[Martin Sandve Alnæs]

It's not to avoid evaluating the conditional, it's because otherwise differentiation can lead to (condition ? trial_function: 0) which is not well defined because 0 loses information about the test space. Is (condition ? trial_function: 0) * test_function a linear or bilinear form? It's not well defined.

This could possibly be remedied by first solving this issue:

https://bitbucket.org/fenics-project/ufl/issues/13/avoid-simplifying-0-v-dx-0-dx-by

An alternative that others have used before is to reformulate your conditionals with a Heaviside function.

[Jan Blechta]

Is (condition ? trial_function: 0) * test_function a linear or bilinear form?

From mathematical point of view it is bilinear form if condition does not depend on trial_function.

[Martin Sandve Alnæs]

Poor formulation from my side, I didn't mean to suggest otherwise. It's a software problem.

[Tanyakarn Treeratanaphitak]

I can't seem to find documentation on implementing Heaviside functions in FEniCS. Would you be able to point me in the right direction?

Thank you.

[Martin Sandve Alnæs]

Actually ufl/ffc issue, partially fixed and I'll finish the job.

[Martin Sandve Alnæs]

Resolved with latest ufl and ffc.

[Anders Logg (Chalmers)]

Removing milestone: 1.7 (automated comment)