The correct way to compute errors

In this tutorial, we learn that if we are interested in finding the error of a numerical solution u_h to a given problem with exact solution u_{ex}, we should not compute directly the integral

\displaystyle \epsilon := \int_\Omega | u_h - u_{ex} |^2

which corresponds to

V = FunctionSpace(mesh, ("Lagrange", degree))
error = form((uh - u_ex)**2 * ufl.dx)
E = np.sqrt(comm.allreduce(assemble_scalar(error), MPI.SUM))

Question 1: I assume, that form() is the function that calls the quadrature function in the background, which finds an approximation (depending on the quadrature rule and the degree of u_h and u_{ex}) to the integral above, right?

The reason is that after expanding the square root in the integral, the compute tries to integrate

\displaystyle u_{ex}^2 + u_h^2 - 2u_{ex}u_h

where the terms u_{ex}^2 + u_h^2 and 2u_{ex}u_h are of a similar size, if u_{ex} and u_h are of similar size (which is the case when the error is small) - leading to round-off errors.

Instead of that we should interpolate both u_h and u_{ex} into a higher order space of polynomials.

W = FunctionSpace(mesh, (family, degree + degree_raise))
u_W = Function(W)
u_ex_W = Function(W)

and then we subtract the nodal values (the expansion coefficients of the FEM functions) of the two functions and compute the error as before:

e_W = Function(W)
e_W.x.array[:] = u_W.x.array - u_ex_W.x.array
error = form(ufl.inner(e_W, e_W) * ufl.dx)
E = np.sqrt(comm.allreduce(assemble_scalar(error), MPI.SUM))

Question 2: why do we subtract the arrays of nodal values instead of subtracting the Function objects inside the ufl.inner() like we did for the naive error computation?

Question 3: Why should this procedure give a more accurate error? Morally speaking, we are loosing information by interpolating the exact solution in a polynomial space, right? I tried to work it out, but I can not see why this approach avoids round-off errors.

Yes. You can also define the quadrature rule using the metadata option in ufl.dx.

Again, then it leads to the expansion of the inner product.

You avoid round-off errors, at the cost of using the interpolated solution (and not the exact solution at quadrature points).

Ahhh so the main idea is to subtract nodal values, and the interpolation into the high-order space is just a necessary extra step we do to compensate the loss we get from the interpolation!
Thanks a lot, I appreciate your fast response.

How do I compute the error without interpolating the exact solution?

Say I have an exact solution defined like this:

# exact solution
class u_exact:
    def __init__(self):
        self.t = 0.0
    def eval(self, x):
        return np.sin(np.pi*x[0]) + 10*self.t

Can I avoid the interpolation in the high order polynomial space? Can I compute the error using the exact representation?

The reason I am asking is the Gibbs oscillation phenomenon.

If the exact solution can be represented in ufl, you can write it as (taking your example above):

x = ufl.SpatialCoordinate(mesh)
t = dolfinx.fem.Constant(mesh, 0.1)
u_exact = ufl.sin(ufl.pi*x[0]) + 10 * t_c

and use it in your error computation. Then x[0] is evaluated at the quadrature points in the error assembly, and has not been put in a specific function space.