Confirming functional superconvergence

On 10 August 2012 20:06, Jason Hicken
<email address hidden> wrote:
> New question #205519 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/205519
>
> Hi,
>
> First: remarkable piece of software engineering.
>

Thanks.

> Now my problem. I am trying to confirm functional superconvergence where the solution comes from a Galerkin approximation to a simple 1-D Poisson problem. For a given polynomial order p, it is my understanding that an integral functional should be 2p order accurate for this problem. I obtain this rate for p = 1 and p = 2, but not for higher-order polynomials.
>
> The function below can be used to reproduce the effect: e.g., run with (num=20, p = 3) and (num=40, p=3) and then estimate the convergence rate using log(Error_Functional(num=40)/Error_Functional(num=20))/log(0.5).
>
> Any insights would be appreciated.
>

The errornorm function interpolates finite element functions in an FE
space and defaults to degree 3. You can force a higher degree, e.g.

     print "L2 solution error:", errornorm(u0, u, degree = p+6)

I think this explains why you run into an issue for p >2.

Also, DOLFIN/FFC interpolates all functions in a finite element space,
so when computing errors some care is necessary to be sure that
analytical expressions are evaluated with sufficient accuracy. Doing

    u0 = Expression("1.0 + x[0]*sin(pi*x[0])")
    f = Expression("2.0*pi*cos(pi*x[0]) - pi*pi*x[0]*sin(pi*x[0])")

the form compiler FFC will try to determine a suitable order Lagrange
basis. This should be high enough, but to be sure you can force the
degree, say

    u0 = Expression("1.0 + x[0]*sin(pi*x[0])", degree=p+6)
    f = Expression("2.0*pi*cos(pi*x[0]) - pi*pi*x[0]*sin(pi*x[0])",
degree=p+6)

Garth

> Jason
> _____________________________________________________________
>
> from dolfin import *
> import numpy
>
> def solve(num, p):
> """Solves a 1-D BVP and returns the L2 error, the functional and its error
>
> *Inputs*
> num: number of finite-elements
> p: polynomial order
>
> *Output*
> functional error
> functional error estimate
> correct functional error
>
> """
>
> # Create mesh and define function space
> mesh = UnitInterval(num)
> V = FunctionSpace(mesh, "CG", p)
>
> # Define Dirichlet boundary (x = 0 or x = 1)
> def boundary(x):
> return x[0] < DOLFIN_EPS or x[0] > 1.0 - DOLFIN_EPS
> def u0_boundary(x, on_boundary):
> return on_boundary
>
> # Define exact solution, source, and boundary conditions
> u0 = Expression("1.0 + x[0]*sin(pi*x[0])")
> f = Expression("2.0*pi*cos(pi*x[0]) - pi*pi*x[0]*sin(pi*x[0])")
> bc = DirichletBC(V, u0, u0_boundary)
>
> # define the bilinear form a and linear form L
> v = TestFunction(V)
> u = TrialFunction(V)
> a = -inner(grad(u), grad(v))*dx
> L = f*v*dx
>
> # Solve variational problem and print L2 solution error
> u = Function(V)
> problem = LinearVariationalProblem(a, L, u, bcs=bc)
> solver = LinearVariationalSolver(problem)
> solver.solve()
> print "L2 solution error:", errornorm(u0, u)
>
> # Compute functional and its error
> energy = u*dx
> E = assemble(energy)
> E_e = 1.0/pi + 1.0
> print "Functional value:", E
> print "Functional error:", E_e - E
>
> return [E, E_e - E]
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.