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!
