Malfunctioning of errornorm?

Question

Hello.

I noticed that the errornorm function does not work well in my test code.
In the error computation part, I computed H1 errors between u and u_exact in two different ways, one with the errornorm function and the other with manual formula. It looks that the error obtained by errornorm is not correct.

For your convenience, I attached its result here.

Error norm of u_H1_error
mesh size =0.2500000 Ei=0.2709354 Ei+1=0.2789436 r=-0.04
mesh size =0.1250000 Ei=0.2789436 Ei+1=0.3101136 r=-0.15
mesh size =0.0625000 Ei=0.3101136 Ei+1=0.3196966 r=-0.04
mesh size =0.0312500 Ei=0.3196966 Ei+1=0.3232146 r=-0.02

Error norm of u_H1man_error
mesh size =0.2500000 Ei=0.0623376 Ei+1=0.0140184 r=2.15
mesh size =0.1250000 Ei=0.0140184 Ei+1=0.0032957 r=2.09
mesh size =0.0625000 Ei=0.0032957 Ei+1=0.0007966 r=2.05
mesh size =0.0312500 Ei=0.0007966 Ei+1=0.0001957 r=2.03

Is there anyone who know the reasons?

from dolfin import *
import numpy

h = []  # mesh sizes
E = []  # errors

f=Expression('sin(pi*x[1])*(pi*pi*x[0])', degree = 10)

for nx in [4,8,16,32,64]:
    h.append(1.0/nx)
    mesh = UnitSquareMesh(nx, nx)
    Vc = FiniteElement("CG", mesh.ufl_cell(), 2)
    V0 = FiniteElement("DG", mesh.ufl_cell(), 0)
    Vh = FunctionSpace(mesh, MixedElement([Vc, V0]))

    u_exact=Expression('x[0]*sin(pi*x[1])', degree = 6, domain = mesh)

    u, u0 = TrialFunction(Vh)
    v, v0 = TestFunction(Vh)
    n = FacetNormal(mesh)
    hcell = CellSize(mesh)
    a = inner(grad(u), grad(v))*dx \
    + (- inner(avg(grad(u)), n("+"))*jump(v0) - inner(avg(grad(v)), n("+"))*jump(u0) \
    +  (10.0/avg(hcell))*jump(u0)*jump(v0) )*dS \
    + (- inner(grad(u), n)*(v+v0) - inner(grad(v), n)*(u+u0))*ds + (10.0/hcell)*(u+u0)*(v+v0)*ds
    b = inner(f, v)*dx + inner(f, v0)*dx \
    + 10./hcell*u_exact*(v+v0)*ds - u_exact*dot(grad(v), n)*ds

    u_sol = Function(Vh)
    solve(a == b, u_sol)
    (u, u0) = u_sol.split(deepcopy=True)

    # Compute errors        
    uH1_error= errornorm(u_exact, u, norm_type = 'H1', degree_rise=3)
    uH1man_error= sqrt(assemble(dot(grad(u_exact) - grad(u), grad(u_exact) - grad(u))*dx))
    errors = {'u_H1_error' : uH1_error, 'u_H1man_error' : uH1man_error}
    E.append(errors)

# Print convergence rates
from math import log as ln
error_types = E[0].keys()
for error_type in sorted(error_types):
    j = len(E)
    print '\nError norm of', error_type
    for i in range(0, j-1):
        r = ln(E[i+1][error_type]/E[i][error_type])/ln(0.5)  # E is a list of errors
        print 'mesh size =%.7f Ei=%.7f Ei+1=%.7f r=%.2f' % (h[i], E[i][error_type], E[i+1][error_type], r)

MiroK · Answer 1 · Jan 13, 2017

Hi John, your manual norm is really $H^1_0$, right? With

uH1_error= errornorm(u_exact, u, norm_type = 'H10', degree_rise=3)

both norms give second order convergence,

answered Jan 13, 2017 by MiroK FEniCS Expert (80,920 points)

Hi, Miro. I computed only the gradient part of u without u itself. I don't see why the manual norm is H_0^1.

commented Jan 13, 2017 by JeonghunLee FEniCS User (1,050 points)

Sorry about the confusion, above I meant $H^1$ seminorm. This is actually what erronorm
with norm_type='H10' computes, see here.

commented Jan 13, 2017 by MiroK FEniCS Expert (80,920 points)

I see. The two norms in the code are different. Thanks!

commented Jan 13, 2017 by JeonghunLee FEniCS User (1,050 points)

Malfunctioning of errornorm?

1 Answer