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)