Using conditionals to prevent divide by zero in residual

Question

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.

My expression is:

$$ u = 0 \quad ||u|| \leq 0.001$$
$$ u = (2/||u||)(1 + ||u||^{0.7}) \quad 10 > ||u|| > 0.001$$
$$ u = 0.14 \quad ||u|| \geq 10 $$

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 - true/false) * (expression) and since I have something divided by zero, it does not matter what the first part returns, the result will still be NaN. Is this what is happening? If yes, is there a way around this that does not involve modifying the header file generated by ffc? This does not seem to be an issue when I tried using the same conditional statement with the linear solver, but I need the nonlinear solver for my application.

Thank you.

Answer 1

would a simple if-statement be ok for you?

In that case you could try something like

l2norm = assemble(inner(u,u)) # computes actual value of the norm.
if l2norm < tolerance:
    f = Constant(0.0)
elseif :
    f = sime_other_form

regards

Comments

The expression in question is in the variational problem, I can't use if statements there.

---

Hi threerat,
I am facing the same problem. Did you already manage to solve it?
Thanks in advance,
Santiago

---

The issue seems to be fixed in version 1.6 (see here). I haven't had the chance to test it out though.

Answer 2

I have faced this issue and solved it by supplying a non-zero initial guess as input to the solver.

Answer 3

Another solution is to define

norm = sqrt(dot(u,u) + DOLFIN_EPS)

This will not only avoid division by zero issues, but will also help speeding up the convergence of the Newton method by increasing the regularity of the nonlinear forcing term.