Finding error in computing 'L2' norm in SIPG method

General
1

As I am using SIPG method in -delta u = f for u = 1+x^2+y^2 on unit square mesh
The code corresponding this problem is given as

from dolfin import *
set_log_level(LogLevel.ERROR)

for j in range(1,6):
    nx = 2**(j+1)
    mesh = UnitSquareMesh(nx, nx)
    V = FunctionSpace(mesh, 'DG', 1)
    u = TrialFunction(V)
    v = TestFunction(V)
    n = FacetNormal(mesh)
    h = CellDiameter(mesh)
    h_avg = (h('+') + h('-'))/2
    alpha = 200

    #Define boundary
    u_D = Expression('1+x[0]*x[0]+x[1]*x[1]', degree=4)
    def boundary(x, on_boundary):
        return on_boundary
    bc = DirichletBC(V, u_D, boundary)
    g = Expression('1+x[0]*x[0]+x[1]*x[1]', degree=4) # u = g on boundary = Dirichlet boundary


    #functions
    
    f = Constant('-4')
    a = dot(grad(u),grad(v))*dx - dot(avg(grad(u)),jump(v,n))*dS - dot(avg(grad(v)),jump(u,n))*dS + (alpha/h_avg) * dot(jump(u,n),jump(v,n))*dS - dot(grad(u),v*n)*ds - dot(grad(v),u*n)*ds + (alpha/h) *u*v*ds
    L = f*v*dx + (alpha/h)*g*v*ds
    u = Function(V)
    solve(a==L, u)
    error_L2 = errornorm(u_D, u, 'L2')
    error_H1 = errornorm(u_D, u, 'H1')
    print ('nx=%.2f error_L2=%8.4E error_H1=%8.4E' % ( nx, error_L2, error_H1))

Here as I am increasing mesh refinements, error is increasing .
Can anyone let me know what is the problem in defining a?

2

You are missing the term corresponding to:

i.e.
L = f*v*dx + (alpha/h)*g*v*ds - dot(grad(v), g*n)*ds
which gives you the following errors:

nx=4.00 error_L2=1.2074E-02 error_H1=2.0286E-01
nx=8.00 error_L2=3.0379E-03 error_H1=1.0135E-01
nx=16.00 error_L2=7.6216E-04 error_H1=5.0673E-02
nx=32.00 error_L2=1.9089E-04 error_H1=2.5338E-02
nx=64.00 error_L2=4.7768E-05 error_H1=1.2669E-02
1 Like