Results depend on type of mesh in simple heat equation example

Hi all, I am a new FEniCS user. I am trying to practice the configuration of FEniCS from the simplest one, one of which is a heat equation. I simply copied and pasted the example code, and changed the input mesh. Then, the results were very weird as follows. When I used the crossed mesh, the error has regionally big values whereas the error was small with the left and right mesh. The code is attached here.

FEniCS tutorial demo program: Heat equation with Dirichlet conditions.
Test problem is chosen to give an exact solution at all nodes of the mesh.

u’= Laplace(u) + f in the unit square
u = u_D on the boundary
u = u_0 at t = 0

u = 1 + x^2 + alphay^2 + \betat
f = beta - 2 - 2*alpha

from __future__ import print_function
from fenics import *
import numpy as np

T = 2.0            # final time
num_steps = 10     # number of time steps
dt = T / num_steps # time step size
alpha = 3          # parameter alpha
beta = 1.2         # parameter beta

# Create mesh and define function space
nx = ny = 8
#mesh = UnitSquareMesh(nx, ny,'crossed')
#mesh = UnitSquareMesh(nx, ny,'left')
#mesh = UnitSquareMesh(nx, ny,'right')

V = FunctionSpace(mesh, 'P', 1)

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

def boundary(x, on_boundary):
    return on_boundary

bc = DirichletBC(V, u_D, boundary)

# Define initial value
u_n = interpolate(u_D, V)
#u_n = project(u_D, V)

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(beta - 2 - 2*alpha)

F = u*v*dx + dt*dot(grad(u), grad(v))*dx - (u_n + dt*f)*v*dx
a, L = lhs(F), rhs(F)

# Time-stepping
u = Function(V)
t = 0
vtkfile = File('tutorial.pvd')
for n in range(num_steps):

    # Update current time
    t += dt
    u_D.t = t

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

    # Plot solution
    vtkfile << (u,t)

    # Compute error at vertices
    u_e = interpolate(u_D, V)
    error = np.abs(u_e.vector()[:] - u.vector()[:]).max()
    print('t = %.2f: error = %.3g' % (t, error))

    # Update previous solution

# Hold plot

A specific description of code is given here,

This is the error plot over the domain with the crossed mesh

What you are observing is a super-convergence for structured grids.
The “left” and “right” orientation gives the approximation field certain nice properties, which a “crossed” mesh does not.

You can observe this by for instance increasing the polynomial degree of your approximation space from 1 to 2. Then you will see a maximum error of order 1e-14, which is close to machine precision.
You should also consider the L^2(\Omega) error as opposed to the max error at vertices, i.e.

print("t", errornorm(u_D, u))

Which gives you a better picture of the accumulated error.

Thank you for the considerate comment. But I have never heard about super-convergence in this field. Can you elaborate a little bit more about it please??

The meshes with “left”/“right” orientation, combined with a P1 finite element perfectly approximates the solution to your quadratic problem at the mesh vertices.

However, as you are using a first order function space to approximate a second order polynomial, the solution is not exact at any other point. That is why one usually uses the “L2” error computation, see Solving PDEs in Minutes - <br> The FEniCS Tutorial Volume I and the remaining bits of the chapter.

You should also consider computing convergence rates, as described in: Solving PDEs in Minutes - <br> The FEniCS Tutorial Volume I