Results depend on type of mesh in simple heat equation example

1

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
    plot(u)
    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
    u_n.assign(u)

# Hold plot
#interactive()

A specific description of code is given here,

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

2

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.

3

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??

4

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

2 Likes