Got nan value when interpolate an expression

dolfinx
,
1

I found that expression \nabla|\nabla f| will generate nan value in some cases. Below is an example. As I’m new to the discourse site, I cannot upload a picture. But the output in the code shows nan value in paraview.

Dolfinx version: 0.10.0

from pathlib import Path

from mpi4py import MPI

import ufl
from basix.ufl import element, mixed_element
from dolfinx import default_real_type
from dolfinx.io import VTXWriter
from dolfinx.fem import Function, Expression, functionspace
from dolfinx.mesh import create_rectangle

from ufl import grad, sqrt, dot

msh = create_rectangle(MPI.COMM_WORLD, ((-1, -1), (1, 1)), [500, 500])
P1 = element("Lagrange", msh.basix_cell(), 2, dtype=default_real_type)
ME = functionspace(msh, P1)

u0 = Function(ME)
u1 = Function(ME)

u0.x.array[:] = 0.0
u1.interpolate(Expression(grad(sqrt(dot(grad(u0), grad(u0))))[0], ME.element.interpolation_points))
u1.x.scatter_forward()

folder = Path("test")
bpfile = VTXWriter(MPI.COMM_WORLD, folder / "data.bp", u1)
bpfile.write(0.0)
bpfile.close()

While this problem can be fixed by using an intermediate function:

from pathlib import Path

from mpi4py import MPI

import ufl
from basix.ufl import element, mixed_element
from dolfinx import default_real_type
from dolfinx.io import VTXWriter
from dolfinx.fem import Function, Expression, functionspace
from dolfinx.mesh import create_rectangle

from ufl import grad, sqrt, dot

msh = create_rectangle(MPI.COMM_WORLD, ((-1, -1), (1, 1)), [500, 500])
P1 = element("Lagrange", msh.basix_cell(), 2, dtype=default_real_type)
ME = functionspace(msh, P1)

u0 = Function(ME)
u1 = Function(ME)
f = Function(ME)

u0.x.array[:] = 0.0
f.interpolate(Expression(sqrt(dot(grad(u0), grad(u0))), ME.element.interpolation_points))
f.x.scatter_forward()
u1.interpolate(Expression(grad(f)[0], ME.element.interpolation_points))
u1.x.scatter_forward()

folder = Path("test")
bpfile = VTXWriter(MPI.COMM_WORLD, folder / "data.bp", u1)
bpfile.write(0.0)
bpfile.close()```
2

When you take the gradient of this expression, you end up with something that is infinity in the denominator, i.e.
\nabla\left( \sqrt{g(x)} \right)= \frac{1}{\sqrt{g(x)}}g'(x), and in your case g(x)=0 thus this goes to infinity.
To avoid this singularity you can make the sqrt conditional:

g = dot(grad(u0), grad(u0))
cond_expr = g > 0
g_cond = conditional(cond_expr, g, 1e-32)
expr = grad(sqrt(g_cond))[0]