Initial thought: could you be generating a NaN in your own code as the mesh is refined?  E.g. if you're interpolating some expression that leads to a 0/0 at some points in the domain.

[Buesing, Henrik] I am dividing by Delta_h_facet (see [1]). So for the cell in the origin, the node should move closer to zero. But, since I am using DQ0 it should not be zero, right?

I am also dividing by dt. Should I reformulate s.t. I multiply by dt?

Thank you!

Henrik

[1]

x_func_expr = Expression("x[0]")

y_func_expr = Expression("x[1]")

z_func_expr = Expression("x[2]")

x_func = interpolate(x_func_expr, DG1)

y_func = interpolate(y_func_expr, DG1)

z_func = interpolate(z_func_expr, DG1)

Delta_h_facet = sqrt(jump(x_func)**2

+ jump(y_func)**2

                   + jump(z_func)**2)