Poisson solution with fine mesh gets worse and changes on each run

1

jupyter-notebook version 6.0.3
on Ubuntu 20.04
FEniCS version 2019.1.0

The code from the tutorial “Solving the Poisson equation” behaves strangely when I increase the resolution of the mesh.

Firstly, the errors seem to increase (or at least fluctuate) with finer meshes:

# mesh 8x8
error_L2  = 0.008235098073354848
error_max = 1.3322676295501878e-15
# mesh 16x16
error_L2  = 2.1118912621708272
error_max = 5.329070518200751e-15
#mesh 32x32
error_L2  = 2.0750004430294164
error_max = 2.1760371282653068e-14

But the strangest thing happens at a 64x64 mesh: The solution changes on each run and the errors are huge:

# mesh 64
error_L2  = 6.071644593466394
error_max = 26.340951344432604
# mesh 64 again
error_L2  = 3.904042361757662
error_max = 15.981625955363619
# mesh 64 again
error_L2  = 11.473897896657363
error_max = 46.409562760391076

This is also visible in the plot, where we seem to see a changing cluster of random extrema (yellow/blue spots) and incorrectly constant boundary values (green near the boundary).

Aside from the mesh resolution, my code is unchanged from the tutorial page (my actual problem is a cubic polynomial, but the strange behaviour is the same):

from fenics import *

# Create mesh and define function space
mesh = UnitSquareMesh(8, 8)
V = FunctionSpace(mesh, 'P', 1)

# Define boundary condition
u_D = Expression('1 + x[0]*x[0] + 2*x[1]*x[1]', degree=2)

def boundary(x, on_boundary):
    return on_boundary

bc = DirichletBC(V, u_D, boundary)

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(-6.0)
a = dot(grad(u), grad(v))*dx
L = f*v*dx

# Compute solution
u = Function(V)
solve(a == L, u, bc)

# Plot solution and mesh
plot(u)

# Save solution to file in VTK format
vtkfile = File('poisson/solution.pvd')
vtkfile << u

# Compute error in L2 norm
error_L2 = errornorm(u_D, u, 'L2')

# Compute maximum error at vertices
vertex_values_u_D = u_D.compute_vertex_values(mesh)
vertex_values_u = u.compute_vertex_values(mesh)
import numpy as np
error_max = np.max(np.abs(vertex_values_u_D - vertex_values_u))

# Print errors
print('error_L2  =', error_L2)
print('error_max =', error_max)

(A similar error was described here with no solution.)

2

I cannot reproduce your error on my ubuntu machine. With some small modifications to your code:

from fenics import *


l2e = []
for n in range(6):
    # Create mesh and define function space
    mesh = UnitSquareMesh(2**(n+1), 2**(n+1))
    V = FunctionSpace(mesh, 'P', 1)

    # Define boundary condition
    u_D = Expression('1 + x[0]*x[0] + 2*x[1]*x[1]', degree=2)

    def boundary(x, on_boundary):
        return on_boundary

    bc = DirichletBC(V, u_D, boundary)

    # Define variational problem
    u = TrialFunction(V)
    v = TestFunction(V)
    f = Constant(-6.0)
    a = dot(grad(u), grad(v))*dx
    L = f*v*dx

    # Compute solution
    u = Function(V)
    solve(a == L, u, bc)

    # Plot solution and mesh
    plot(u)

    # Save solution to file in VTK format
    vtkfile = File('poisson/solution.pvd')
    vtkfile << u

    # Compute error in L2 norm
    error_L2 = errornorm(u_D, u, 'L2')
    l2e.append(error_L2)

    # Compute maximum error at vertices
    vertex_values_u_D = u_D.compute_vertex_values(mesh)
    vertex_values_u = u.compute_vertex_values(mesh)
    import numpy as np
    error_max = np.max(np.abs(vertex_values_u_D - vertex_values_u))

    # Print errors
    print('error_L2  =', error_L2)
    print('error_max =', error_max)

import numpy as np
l2e = np.array(l2e, dtype=np.double)

print(np.log(l2e[:-1]/l2e[1:])/np.log(2.0))

I see optimal convergence rates [2. 2. 2. 2. 2.] (i.e. a second order method).