Evaluate errornorm at points that are different from the mesh

1

Hi, I am new to FEniCS and trying to solve the bourdary problem of the Poissson Equation in the tutorial (Chapter 2: https://fenicsproject.org/pub/tutorial/html/._ftut1004.html#ch:fundamentals). Note I am using 2019.2.0.dev0 of dolphin and ufl package.

I could reproduce the result, but I could not understand how FEniCS calculates the L2 erorr. Using the code below, I obtained the L2 error of 0.008235098073354943 and max error at mesh points of 1.3322676295501878e-15. I think the max error is correct, while the L2 error is pretty large. My understanding is that the FEniCS interpolates the finite element solution and evaluates the L2 norm analytically over the domain rather than evaluating it on the mesh points used in the finite element solution. Is this correct?

import dolfin as dl
import ufl

import math
import numpy as np
import logging

import matplotlib.pyplot as plt
%matplotlib inline

import sys
import os
sys.path.append( os.environ.get('HIPPYLIB_BASE_DIR', "../") )
from hippylib import nb

logging.getLogger('FFC').setLevel(logging.WARNING)
logging.getLogger('UFL').setLevel(logging.WARNING)
dl.set_log_active(False)

n = 8
mesh = dl.UnitSquareMesh(n , n)
degree = 1
nb.plot(mesh)

Vh = dl.FunctionSpace(mesh, 'Lagrange', degree)
print("dim(Vh) = ", Vh.dim())

u = dl.TrialFunction(Vh)
v = dl.TestFunction(Vh)

tol = dl.DOLFIN_EPS

formula = '1 + x[0]*x[0] + 2 * x[1]*x[1]'
u_D = dl.Expression(formula, degree = 2)
tol = dl.DOLFIN_EPS
def boundary(x):
    return abs(x[0]) < tol or abs(x[1]) < tol or abs(x[0] - 1) < tol or abs(x[1] - 1) < tol
bc = dl.DirichletBC(Vh, u_D, boundary)

f = dl.Expression('-6', degree = 0)
a = ufl.inner(ufl.grad(u), ufl.grad(v))*ufl.dx
L = f * v * ufl.dx

uh = dl.Function(Vh)
dl.solve(a == L, uh, bc)
# A, b = dl.assemble_system(a, L, bcs = bc)
# dl.solve(A, uh.vector(), b, 'cg')

nb.plot(uh)
nb.plot(mesh)

vertex_values_u_D = u_D.compute_vertex_values(mesh)
vertex_values_uh = uh.compute_vertex_values(mesh)
error_max = np.max(np.abs(vertex_values_u_D - vertex_values_uh))

error_L2 = dl.errornorm(u_D, uh, 'L2')

# err_L2 = math.sqrt(dl.assemble((u_D - uh) ** 2*ufl.dx))

print(error_L2)
print(error_max)

Also, I want to evaluate the finite element solution at arbitrary points. Is there any way to get interpolated solution and evaluate it on points I specify? I tried the code below, but I got an error. I appreciate any help!

tol = 0.001
y = np.linspace(-1 + tol, 1 - tol, 101)
points = [(0, y_) for y_ in y]
u_line = np.array([u(point) for point in points])
plt.plot(y, u_line)

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
<ipython-input-117-69983c27e77c> in <module>
      2 y = np.linspace(-1 + tol, 1 - tol, 101)
      3 points = [(0, y_) for y_ in y]
----> 4 u_line = np.array([u(point) for point in points])
      5 plt.plot(y, u_line)

<ipython-input-117-69983c27e77c> in <listcomp>(.0)
      2 y = np.linspace(-1 + tol, 1 - tol, 101)
      3 points = [(0, y_) for y_ in y]
----> 4 u_line = np.array([u(point) for point in points])
      5 plt.plot(y, u_line)

/usr/lib/petsc/lib/python3/dist-packages/dolfin/function/function.py in __call__(self, *args, **kwargs)
    353 
    354         # The actual evaluation
--> 355         self._cpp_object.eval(values, x)
    356 
    357         # If scalar return statement, return scalar value.

RuntimeError: 

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to evaluate function at point.
*** Reason:  The point is not inside the domain. Consider calling "Function::set_allow_extrapolation(true)" on this Function to allow extrapolation.
*** Where:   This error was encountered inside Function.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.2.0.dev0
*** Git changeset:  unknown
*** -------------------------------------------------------------------------
2

Remember that the FE approximation u_h(x) is actually a continuous field. Consider the L_2 norm
\Vert u - u_h \Vert_{L_2(\Omega)} = \sqrt{\int_\Omega (u - u_h)^2 \; \mathrm{d}x},
this gives us information about the error in the FE approximation at all points in our domain. The absolute value of this measure may not be as important as the relative value compared with the same measure of the function. For example, \Vert u - u_h \Vert_{L_2(\Omega)} may be large; however, \Vert u - u_h \Vert_{L_2(\Omega)}/\Vert u\Vert_{L_2(\Omega)} may be small.

Futhermore there are lots of numerical analysis tools we may exploit to learn more about how this error should converge, and at what rates. This equips us with avenues for verification of our FE methods and code.

What you propose by checking the solution “at points” is a discrete \ell_1 vector norm at arbitrary points. This throws away information regarding the accuracy of the FE problem in the entire domain. Furthermore we cannot exploit the numerical toolbox offered by Galerkin orthogonality for further analysis.

Finally the reason you see exact precision at the nodes is because your analytical solution, u = 1+x^2+2y^2 can be exactly represented in a the FE space V^{h,2}. You would’ve learned this by increasing the polynomial degree by one degree=2 and observing that in this case \Vert u - u_h \Vert_{L_2(\Omega)} = 0 to machine precision. You couldn’t have known this by employing the discrete \ell_1 norm at arbitrary points.

Hughes 2012 is a good introduction to the mathematics of finite elements for engineers and scientists.

1 Like
3

Hi Nate,

Thank you very much for your quick and in-depth explanations! As you mentioned, I should have used the relative error to compare the result with other numerical methods.

In terms of the last part, it perfectly makes sense. When I use degree = 2, I got a pretty good solution (around e-14). Also, thank you for the reference!

1 Like
4

To add to @nate’s reply, your second issue is quite easy to explain, as your domain

Goes from [0,1]x[0,1], while your y coordinate goes from [-1,1]

Resulting in the very clear error message

Note that in chapter 2.3.11 of the accompanying FEniCS tutorial book they briefly cover error estimates.

2 Likes
5

Hi dokken,

Wow, why did not I notice it…! Thank you so much!