This is a read only copy of the old FEniCS QA forum. Please visit the new QA forum to ask questions

Using conditionals to prevent divide by zero in residual

0 votes

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.

asked Jul 2, 2015 by ttreerat FEniCS Novice (220 points)

3 Answers

0 votes

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

answered Aug 6, 2015 by multigrid202 FEniCS User (3,780 points)

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.

0 votes

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

answered Feb 3, 2017 by KristianE FEniCS Expert (12,900 points)
0 votes

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.

answered Feb 3, 2017 by umberto FEniCS User (6,440 points)
...