Unexpectedly high L2 error when comparing coarse and fine resolution approximation

Dear community,

I observed an unexpectedly high L2 error when I compared the approximated solution of the poisson problem considered in the error control tutorial with a high-resolution approximation instead of the analytical solution. Unfortunately, I can’t explain why the error is so much higher than when comparing with the analytical solution. Do you have a reason for this?

The reason why I want to estimate the L2 error with respect to a high-resolution approximation is that I do not have an analytical solution for my use case.

Here is the code:

from dolfinx.fem import (
    Expression,
    Function,
    FunctionSpace,
    assemble_scalar,
    dirichletbc,
    form,
    locate_dofs_topological,
)
from dolfinx.fem.petsc import LinearProblem
from dolfinx.mesh import create_unit_square, locate_entities_boundary

from mpi4py import MPI
from ufl import SpatialCoordinate, TestFunction, TrialFunction, div, dx, grad, inner

import ufl
import numpy as np


def u_ex(mod):
    return lambda x: mod.cos(2 * mod.pi * x[0]) * mod.cos(2 * mod.pi * x[1])


u_numpy = u_ex(np)
u_ufl = u_ex(ufl)


def solve_poisson(N=10, degree=1):
    mesh = create_unit_square(MPI.COMM_WORLD, N, N)
    x = SpatialCoordinate(mesh)
    f = -div(grad(u_ufl(x)))
    V = FunctionSpace(mesh, ("Lagrange", degree))
    u = TrialFunction(V)
    v = TestFunction(V)
    a = inner(grad(u), grad(v)) * dx
    L = f * v * dx
    u_bc = Function(V)
    u_bc.interpolate(u_numpy)
    facets = locate_entities_boundary(
        mesh, mesh.topology.dim - 1, lambda x: np.full(x.shape[1], True)
    )
    dofs = locate_dofs_topological(V, mesh.topology.dim - 1, facets)
    bcs = [dirichletbc(u_bc, dofs)]
    default_problem = LinearProblem(
        a, L, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"}
    )
    return default_problem.solve(), u_ufl(x)


def error_L2(uh, u_ex, degree_raise=3):
    # Create higher order function space
    degree = uh.function_space.ufl_element().degree()
    family = uh.function_space.ufl_element().family()
    mesh = uh.function_space.mesh
    W = FunctionSpace(mesh, (family, degree + degree_raise))
    # Interpolate approximate solution
    u_W = Function(W)
    u_W.interpolate(uh)

    # Interpolate exact solution, special handling if exact solution
    # is a ufl expression or a python lambda function
    u_ex_W = Function(W)
    if isinstance(u_ex, ufl.core.expr.Expr):
        u_expr = Expression(u_ex, W.element.interpolation_points())
        u_ex_W.interpolate(u_expr)
    else:
        u_ex_W.interpolate(u_ex)

    # Compute the error in the higher order function space
    e_W = Function(W)
    e_W.x.array[:] = u_W.x.array - u_ex_W.x.array

    # Integrate the error
    error = form(ufl.inner(e_W, e_W) * ufl.dx)
    error_local = assemble_scalar(error)
    error_global = mesh.comm.allreduce(error_local, op=MPI.SUM)
    return np.sqrt(error_global)


num_elements_list = [4, 8, 16]
u_fine, _ = solve_poisson(N=128)

print("Compare to exact solution")
for num_elements in num_elements_list:
    u_approx, _ = solve_poisson(N=num_elements)
    print(f"{num_elements}: {error_L2(u_approx, u_numpy)}")

print("Compare to finely resolved approximation")
for num_elements in num_elements_list:
    u_approx, _ = solve_poisson(N=num_elements)
    print(f"{num_elements}: {error_L2(u_approx, u_fine)}")

print("L2-error of finely resolved approximation")
print(error_L2(u_fine, u_numpy))

This is the corresponding output:

Compare to exact solution
4: 0.24336088241081913
8: 0.07959834507615321
16: 0.021454239435480606
Compare to finely resolved approximation
4: 1.0979738428201873
8: 1.0589157556605937
16: 0.9489944941705492
L2-error of finely resolved approximation
0.0003439221768301646

I ran the code in a docker container built from dolfinx/dolfinx:stable.

Thank you in advance for your help! Greetings!

Could you specify what version of DOLFINx you use for this code? as it seems like you are here relying on interpolation of solutions between meshes (from v0.6.0).
I would also suggest you store the interpolated solution u_W to file to visualize how it looks, as it might give you some points as to why the error is wrong.

First of all, thanks for the quick reply!

I use DOLFINx version 0.7.0. These are the related FEniCSx package versions installed in the docker container:

fenics-basix 0.7.0
fenics-dolfinx 0.7.0
fenics-ffcx 0.7.0
fenics-ufl 2023.2.0

I have visualized both the highly resolved approximation (128 elements per dimension) as well as the exact solution using ParaView. In order to visualize the exact solution, I interpolated the NumPy function to the same function space I used for the highly resolved approximation. As I am only allowed to upload one figure per post, the exact solution will follow in another post.

Unfortunately, I could not apply a filter to visualize the differences of both plots in ParaView. At first glance, however, no relevant deviations can be seen, which also matches the small L2 error of the high-resolution approximation compared to the exact solution.

Highly resolved approximation

This gives me the impression that the high-resolution approximation is not the reason for the high L2 error. What could be other sources for the high L2 error? Many thanks in advance for your advice. Greetings!

Exact solution

You are not handling the interpolation between meshes correctly.

You should inspect u_ex_W, which give you pointers to what is wrong.
Here is a corrected script that uses the correct handling for non-matching mesh interpolation:

from dolfinx.fem import (
    Expression,
    Function,
    FunctionSpace,
    assemble_scalar,
    dirichletbc,
    form,
    locate_dofs_topological,
)
from dolfinx.io import VTXWriter
from dolfinx.fem.petsc import LinearProblem
from dolfinx.mesh import create_unit_square, locate_entities_boundary

from mpi4py import MPI
from ufl import SpatialCoordinate, TestFunction, TrialFunction, div, dx, grad, inner

import ufl
import numpy as np


def u_ex(mod):
    return lambda x: mod.cos(2 * mod.pi * x[0]) * mod.cos(2 * mod.pi * x[1])


u_numpy = u_ex(np)
u_ufl = u_ex(ufl)


def solve_poisson(N=10, degree=1):
    mesh = create_unit_square(MPI.COMM_WORLD, N, N)
    x = SpatialCoordinate(mesh)
    f = -div(grad(u_ufl(x)))
    V = FunctionSpace(mesh, ("Lagrange", degree))
    u = TrialFunction(V)
    v = TestFunction(V)
    a = inner(grad(u), grad(v)) * dx
    L = f * v * dx
    u_bc = Function(V)
    u_bc.interpolate(u_numpy)
    facets = locate_entities_boundary(
        mesh, mesh.topology.dim - 1, lambda x: np.full(x.shape[1], True)
    )
    dofs = locate_dofs_topological(V, mesh.topology.dim - 1, facets)
    bcs = [dirichletbc(u_bc, dofs)]
    default_problem = LinearProblem(
        a, L, bcs=bcs, petsc_options={"ksp_type": "preonly", "pc_type": "lu"}
    )
    return default_problem.solve(), u_ufl(x)


def error_L2(uh, u_ex, degree_raise=3):
    # Create higher order function space
    degree = uh.function_space.ufl_element().degree()
    family = uh.function_space.ufl_element().family()
    mesh = uh.function_space.mesh
    W = FunctionSpace(mesh, (family, degree + degree_raise))
    # Interpolate approximate solution
    u_W = Function(W)
    u_W.interpolate(uh)

    # Interpolate exact solution, special handling if exact solution
    # is a ufl expression or a python lambda function
    u_ex_W = Function(W)

    if isinstance(u_ex, Function):
        if u_ex.function_space.mesh != mesh:
            from dolfinx.fem import create_nonmatching_meshes_interpolation_data
            interpolation_data = create_nonmatching_meshes_interpolation_data(
                W.mesh._cpp_object,
                W.element,
                u_ex.function_space.mesh._cpp_object)
            u_ex_W.interpolate(u_ex, nmm_interpolation_data=interpolation_data)
        else:
            u_ex_W.interpolate(u_ex)
    elif isinstance(u_ex, ufl.core.expr.Expr):
        u_expr = Expression(u_ex, W.element.interpolation_points())
        u_ex_W.interpolate(u_expr)
    else:
        u_ex_W.interpolate(u_ex)

    u_ex_W.name = "u_ex_W"
    u_W.name = "u_W"
    with VTXWriter(mesh.comm, f"test_u_{family}_{degree}.bp", [u_ex_W, u_W], engine="BP4") as vtx:
        vtx.write(0.0)

    # Compute the error in the higher order function space
    e_W = Function(W)
    e_W.x.array[:] = u_W.x.array - u_ex_W.x.array

    # Integrate the error
    error = form(ufl.inner(e_W, e_W) * ufl.dx)
    error_local = assemble_scalar(error)
    error_global = mesh.comm.allreduce(error_local, op=MPI.SUM)
    return np.sqrt(error_global)


num_elements_list = [8, 16, 32]
u_fine, _ = solve_poisson(N=128)

print("Compare to exact solution")
for num_elements in num_elements_list:
    u_approx, _ = solve_poisson(N=num_elements)
    print(f"{num_elements}: {error_L2(u_approx, u_numpy)}")

print("Compare to finely resolved approximation")
for num_elements in num_elements_list:
    u_approx, _ = solve_poisson(N=num_elements)
    print(f"{num_elements}: {error_L2(u_approx, u_fine)}")

print("L2-error of finely resolved approximation")
print(error_L2(u_fine, u_numpy))

returning

Compare to exact solution
8: 0.0795983450761532
16: 0.021454239435480596
32: 0.005469005005656221
Compare to finely resolved approximation
8: 0.07929481752860122
16: 0.02114171637089746
32: 0.005221454749415455
L2-error of finely resolved approximation
0.0003439221768301677

Thanks for the solution!

Hello dr @dokken I was trying to calculate L2 error between fine mesh solution and coarse mesh solution with this code but unfortunately it’s not working in dolfinx 0.8 version. Can you please help or give me more insight into special handling of non matching mesh solution interpolation and then how to calculate error between the two.

As you haven’t posted what code you are running or what error message you are getting i cannot Give you any further guidance.

The code goes like this. It is throwing attribute error while accessing family for the computed solution. Please let me know if there is any other error in the code.

def error_L2(uh, u_ex, degree_raise=3):
    # Create higher order function space
    degree = uh.function_space.ufl_element().degree
    family = uh.function_space.ufl_element().family()
    mesh = uh.function_space.mesh
    W = functionspace(mesh, (family, degree + degree_raise))
    # Interpolate approximate solution
    u_W = Function(W)
    u_W.interpolate(uh)

    # Interpolate exact solution, special handling if exact solution
    # is a ufl expression or a python lambda function
    u_ex_W = Function(W)

    if isinstance(u_ex, Function):
        if u_ex.function_space.mesh != mesh:
            from dolfinx.fem import create_nonmatching_meshes_interpolation_data
            interpolation_data = create_nonmatching_meshes_interpolation_data(
                W.mesh._cpp_object,
                W.element,
                u_ex.function_space.mesh._cpp_object)
            u_ex_W.interpolate(u_ex, nmm_interpolation_data=interpolation_data)
        else:
            u_ex_W.interpolate(u_ex)
    elif isinstance(u_ex, ufl.core.expr.Expr):
        u_expr = Expression(u_ex, W.element.interpolation_points())
        u_ex_W.interpolate(u_expr)
    else:
        u_ex_W.interpolate(u_ex)

    # Compute the error in the higher order function space
    e_W = Function(W)
    e_W.x.array[:] = u_W.x.array - u_ex_W.x.array

    # Integrate the error
    error = form(ufl.inner(e_W, e_W) * ufl.dx)
    error_local = assemble_scalar(error)
    error_global = mesh.comm.allreduce(error_local, op=MPI.SUM)
    return np.sqrt(error_global)


num_elements_list = [8, 16, 32]
u_fine = solve_adv_diff(N=100)


print("Compare to finely resolved approximation")
for num_elements in num_elements_list:
    u_approx = solve_adv_diff(N=num_elements)
    print(f"{num_elements}: {error_L2(u_approx, u_fine)}")


AttributeError                            Traceback (most recent call last)
<ipython-input-23-74051482dee0> in <cell line: 47>()
     47 for num_elements in num_elements_list:
     48     u_approx = solve_adv_diff(N=num_elements)
---> 49     print(f"{num_elements}: {error_L2(u_approx, u_fine)}")

<ipython-input-23-74051482dee0> in error_L2(uh, u_ex, degree_raise)
      2     # Create higher order function space
      3     degree = uh.function_space.ufl_element().degree
----> 4     family = uh.function_space.ufl_element().family()
      5     mesh = uh.function_space.mesh
      6     W = functionspace(mesh, (family, degree + degree_raise))

AttributeError: '_BasixElement' object has no attribute 'family'  ```

Please format the code with 3x`

```python
# add code here

```

Hello, dokken. I tried the codes you posted here, I get the error as follows
update: I solved this problem by adding a “0.0” in the arguments lists in create_nonmatching_meshes_interpolation_data, and I get the same results as you.

python3 test_inter.py 
Compare to exact solution
8: 0.07959834507615322
16: 0.02145423943548061
32: 0.005469005005656221
Compare to finely resolved approximation
Traceback (most recent call last):
  File "/mnt/g/workdir/MKM/NO_mkm/fenics/test_inter.py", line 108, in <module>
    print(f"{num_elements}: {error_L2(u_approx, u_fine)}")
                             ^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/mnt/g/workdir/MKM/NO_mkm/fenics/test_inter.py", line 68, in error_L2
    interpolation_data = create_nonmatching_meshes_interpolation_data(
                         ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
TypeError: create_nonmatching_meshes_interpolation_data(): incompatible function arguments. The following argument types are supported:
    1. (mesh0: dolfinx.cpp.mesh.Mesh_float32, element0: dolfinx.cpp.fem.FiniteElement_float32, mesh1: dolfinx.cpp.mesh.Mesh_float32, padding: float) -> Tuple[List[int], List[int], List[float], List[int]]
    2. (geometry0: dolfinx.cpp.mesh.Geometry_float32, element0: dolfinx.cpp.fem.FiniteElement_float32, mesh1: dolfinx.cpp.mesh.Mesh_float32, cells: numpy.ndarray[numpy.int32], padding: float) -> Tuple[List[int], List[int], List[float], List[int]]
    3. (mesh0: dolfinx.cpp.mesh.Mesh_float64, element0: dolfinx.cpp.fem.FiniteElement_float64, mesh1: dolfinx.cpp.mesh.Mesh_float64, padding: float) -> Tuple[List[int], List[int], List[float], List[int]]
    4. (geometry0: dolfinx.cpp.mesh.Geometry_float64, element0: dolfinx.cpp.fem.FiniteElement_float64, mesh1: dolfinx.cpp.mesh.Mesh_float64, cells: numpy.ndarray[numpy.int32], padding: float) -> Tuple[List[int], List[int], List[float], List[int]]

Invoked with: <dolfinx.cpp.mesh.Mesh_float64 object at 0x7fe924e87d70>, <dolfinx.cpp.fem.FiniteElement_float64 object at 0x7fe924df01b0>, <dolfinx.cpp.mesh.Mesh_float64 object at 0x7fe924e84410>

Here is the version of fenicsx I installed

fenics-basix              0.7.0           py312h5b9907d_1    conda-forge
fenics-basix-pybind11-abi 0.4.12               h5b9907d_1    conda-forge
fenics-dolfinx            0.7.2           py312had67bdd_103    conda-forge
fenics-ffcx               0.7.0              pyh4af843d_0    conda-forge
fenics-libbasix           0.7.0                hfdc072b_1    conda-forge
fenics-libdolfinx         0.7.2              h962240e_103    conda-forge
fenics-ufcx               0.7.0                h4af843d_0    conda-forge
fenics-ufl                2023.2.0           pyhd8ed1ab_0    conda-forge

Please see the code in updated format.

As @Smith_Jack said, you need to add a float to

I would probably not use 0.0, but 1e-6, see

I am unable to follow his suggestion. The error I am getting comes from attribute error while accessing family for the computed solution. Also, I don’t understand where I need to add float. Can you please edit my code. I will be grateful for your kind response.

Maybe the error results from the inconsistent version of dolfinx. I am using 0.7.2.

I think there is some update for accessing family of a function in 0.8 version, I am using google colab and automatically installs latest update from the website. I am unable to figure out how to access it. Also, Can you please higlight where I need to add float to create_nonmatching_meshes_interpolation_data? is it after the last arguement in my code??

I installed Fenicsx by anaconda, you can change the channel of anaconda to see if there is a appropriate version of fenics. Also, I think the version can be set when you use anaconda to install a package.

As for the argument, it’s the padding argument of the create… function (yes, the last argument), as you can see from the error messages given by python.

1 Like

Thank you for your reply. All of my work is in google colab.Please let me know if you can offer help installing fenicsx in anaconda I am not comfortable with that stuff. Or I will be waiting for dr @dokken’s reply for accessing family for the computed solution in 0.8 version.

How to choose the value for padding in terms of speed and accuracy of interpolation? You gave 1e-6 while the unit test script gave 1e-14, and both of them are correct. However, large padding results in a very slow performance or wrong answer.

This depends on if you run dolfinx with a mesh of float64 or float32 geometry.

In general, running double precision (float64) then a tolerance of 1e-6 would be considered sufficient. The tolerance is needed as the collision detection algorithm we use (GJK) requires squared numbers (norms) to compute distances, meaning that the precision is sqrt(floating_point_epsilon).