Computing Error Norms (C++)

Asked by Pietro Maximoff on 2009-12-29

Hi all

I am trying to compute the error given a known solution. (C++). This is what I did.

I created two .ufl files; the first to genrate a functional:

Error.ufl
=======
element = FiniteElement("Lagrange", "triangle", 3)
v = Coefficient(element)
w = Coefficient(element)
M= ((w - v) * (w - v)) * dx

and the second to represent a higher order function space onto which the exact solution will be interpolated:

HigherOrder.ufl
=============
element = FiniteElement("Lagrange", "triangle", 3)

Then in main.cpp, I had

Function u_exact_high(V_high_order);
u_exact_high.interpolate(u_exact);
Error::functional M(mesh, u_exact_high, u);
double error = std::sqrt(assemble(M));

1. is this the correct way to do it?

2. Is there a function or class that I can use? [I've searched but couldn't find one].

3. Do I follow the same procedure to calculate the H1 seminorm and other norms?

Best wishes

Pietro

Question information

Language:
English Edit question
Status:
Solved
For:
DOLFIN Edit question
Assignee:
No assignee Edit question
Solved by:
Andy R Terrel
Solved:
2009-12-30
Last query:
2009-12-30
Last reply:
2009-12-30
Phil Marinier (lonewolf-13p) said : #1

This is the way I did it. It worked for me. I use something different now, because we decided to calculate the error at a different point in our process, but when I was using this it worked fine.

element = FiniteElement("Lagrange", triangle, 1)

uh = Function(element)
u = Function(element)
e = uh - u

M = dot(e, e)*dx

then in the code:

 ENorm::Functional M;

 M.u = Phi0;
 M.uh = Phi1;

 double errorSqr = assemble(M);
 double error = sqrt(errorSqr);

where Phi0 and Phi1 are dolfin::Functions that I had previously solved for. They are both defined over the same function space.

This is obviously least squares. There are built in methods to find error norms if you use python, but as far as I know, the c++ interface hasn't been built yet. However, if you can find the code in the dolfin files, you can use parts of it. That's what I was doing when I was experimenting with different error norms.

I would also like to know if there is a c++ built in way to do this.

I did this a while ago, not sure if the update has any effect on the code.

This isn't a full answer but hopefully it helps.

Best Andy R Terrel (andy-terrel) said : #2

FFC will take care of interpolation of the computed function for you. So you combine your two forms into one:

Error.ufl
=======
e3 = FiniteElement("Lagrange", "triangle", 3)
e10 = FiniteElement("Quadrature", "triangle", 10)
v = Coefficient(e3)
w = Coefficient(e10)
M= ((w - v) * (w - v)) * dx

This can be done for any functional (so all your seminorms and what not). I don't think this is in dolfin, I think a long time ago there were things like this but because the need to support so many different combinations of finite elements it got yanked.

Pietro Maximoff (segment-x) said : #3

Thanks Andy R Terrel, that solved my question.

Pietro Maximoff (segment-x) said : #4

Also, thanks, Phil, your way was more compact than mine.

Anders Logg (logg) said : #5

On Wed, Dec 30, 2009 at 02:00:23PM -0000, Andy R Terrel wrote:
> Question #95601 on DOLFIN changed:
> https://answers.launchpad.net/dolfin/+question/95601
>
> Status: Open => Answered
>
> Andy R Terrel proposed the following answer:
> FFC will take care of interpolation of the computed function for you.
> So you combine your two forms into one:
>
> Error.ufl
> =======
> e3 = FiniteElement("Lagrange", "triangle", 3)
> e10 = FiniteElement("Quadrature", "triangle", 10)
> v = Coefficient(e3)
> w = Coefficient(e10)
> M= ((w - v) * (w - v)) * dx
>
>
> This can be done for any functional (so all your seminorms and what
> not). I don't think this is in dolfin, I think a long time ago
> there were things like this but because the need to support so many
> different combinations of finite elements it got yanked.

Yes, this is something that is not handled automatically in C++
since it would require pregeneration of large amounts of code to cover
different function spaces.

Andy's solution is good, but you might run into stability problems
(round-off errors) for higher order function spaces.

Look at errornorm.py for a better way to compute the error:

  """Compute the error e = u - uh in the given norm. The parameter k
  denotes the degree of accuracy (degree of piecewise polynomials
  approximating u an uh).

  In simple cases, one may just define

    e = u - uh

  and evalute for example the square of the error in the L2 norm by

    e = u - uh
    assemble(e*e*dx, mesh)

  However, this is not stable w.r.t. round-off errors considering
  that the form compiler may expand the expression above to

    u*u*dx + uh*uh*dx - 2*u*uh*dx

  and this might get further expanded into thousands of terms for
  higher order elements. Thus, the error will be evaluated by adding
  a large number of terms which should sum up to something close to
  zero (if the error is small)."""

--
Anders