DOLFIN

Python function errornorm will give misleading results for higher-order elements

Bug #1036135 reported by Garth Wells
8
This bug affects 1 person
Affects Status Importance Assigned to Milestone
DOLFIN
Fix Released
Undecided
Garth Wells

Bug Description

The Python function 'errornorm' in fem.norms uses a default interpolation order of 3, which will produce wrong results for higher-order problems.

The function should choose an order based on the order of the argument functions, and if this is not possible the user should provide the order.

Revision history for this message
Martin Sandve Alnæs (martinal) wrote :

In the email prior to this report, you say that:
"""
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)
"""

AFAIK this is _not_ correct. When making an Expression("..."), FFC does not even know what the expression string is, and the degree that DOLFIN assigns is independent of the expression string.

If, on the other hand, you use an UFL expression:
    u0 = 1.0 + x[0]*sin(pi*x[0])
and compile a functional
    error = (u-u0)**2*dx
then no interpolation of u0 takes place, since the exact expression for u0 will become part of the code generated by FFC, and the complexity of u0 will be visible to FFC/UFL when estimating an integration order.

Revision history for this message
Garth Wells (garth-wells) wrote : Re: [Bug 1036135] Re: Python function errornorm will give misleading results for higher-order elements

On 13 August 2012 12:35, Martin Sandve Alnæs <email address hidden> wrote:
> In the email prior to this report, you say that:
> """
> 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)
> """
>
> AFAIK this is _not_ correct. When making an Expression("..."), FFC does
> not even know what the expression string is, and the degree that DOLFIN
> assigns is independent of the expression string.
>

Yes, it is independent of the string but DOLFIN does not determine
the degree. FFC makes a guess at the degree. The code is in
ffc/ffc/analysis.py in the function

    _auto_select_degree(elements)

> If, on the other hand, you use an UFL expression:
> u0 = 1.0 + x[0]*sin(pi*x[0])
> and compile a functional
> error = (u-u0)**2*dx
> then no interpolation of u0 takes place, since the exact expression for u0 will become part of the code generated by FFC, and the complexity of u0 will be visible to FFC/UFL when estimating an integration order.
>

With the above, the accuracy issue is shifted from interpolation to
integration. The functional above will be integrated approximately,
and one needs to make sure that the integration is sufficiently
accurate. If u0 is an Expression, the integration is exact but there
is an interpolation error.

Garth

> --
> You received this bug notification because you are a member of DOLFIN
> Core Team, which is subscribed to DOLFIN.
> https://bugs.launchpad.net/bugs/1036135
>
> Title:
> Python function errornorm will give misleading results for higher-
> order elements
>
> To manage notifications about this bug go to:
> https://bugs.launchpad.net/dolfin/+bug/1036135/+subscriptions

Revision history for this message
Anders Logg (logg) wrote :

On Mon, Aug 13, 2012 at 12:27:23PM -0000, Garth Wells wrote:
> On 13 August 2012 12:35, Martin Sandve Alnæs <email address hidden> wrote:
> > In the email prior to this report, you say that:
> > """
> > 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)
> > """
> >
> > AFAIK this is _not_ correct. When making an Expression("..."), FFC does
> > not even know what the expression string is, and the degree that DOLFIN
> > assigns is independent of the expression string.
> >
>
> Yes, it is independent of the string but DOLFIN does not determine
> the degree. FFC makes a guess at the degree. The code is in
> ffc/ffc/analysis.py in the function
>
> _auto_select_degree(elements)
>
>
> > If, on the other hand, you use an UFL expression:
> > u0 = 1.0 + x[0]*sin(pi*x[0])
> > and compile a functional
> > error = (u-u0)**2*dx
> > then no interpolation of u0 takes place, since the exact expression for u0 will become part of the code generated by FFC, and the complexity of u0 will be visible to FFC/UFL when estimating an integration order.
> >
>
> With the above, the accuracy issue is shifted from interpolation to
> integration. The functional above will be integrated approximately,
> and one needs to make sure that the integration is sufficiently
> accurate. If u0 is an Expression, the integration is exact but there
> is an interpolation error.

Another issue here is how the integral is evaluated. The important
point is to not do the following:

  e = \int (u - u_h)^2 = \int u^2 + \int u_h^2 - \int 2*u*uh

since this will lead to a subtraction of two numbers of approximately
the same size (if the error is small).

That is why errornorm interpolates u and u_h into the same space, then
subtracts the dof values, and then performs the integration.

--
Anders

Revision history for this message
Martin Sandve Alnæs (martinal) wrote :

On 13 August 2012 17:35, Anders Logg <email address hidden> wrote:
> On Mon, Aug 13, 2012 at 12:27:23PM -0000, Garth Wells wrote:
>> On 13 August 2012 12:35, Martin Sandve Alnæs <email address hidden> wrote:
>> > In the email prior to this report, you say that:
>> > """
>> > 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)
>> > """
>> >
>> > AFAIK this is _not_ correct. When making an Expression("..."), FFC does
>> > not even know what the expression string is, and the degree that DOLFIN
>> > assigns is independent of the expression string.
>> >
>>
>> Yes, it is independent of the string but DOLFIN does not determine
>> the degree. FFC makes a guess at the degree. The code is in
>> ffc/ffc/analysis.py in the function
>>
>> _auto_select_degree(elements)
>>
>>
>> > If, on the other hand, you use an UFL expression:
>> > u0 = 1.0 + x[0]*sin(pi*x[0])
>> > and compile a functional
>> > error = (u-u0)**2*dx
>> > then no interpolation of u0 takes place, since the exact expression for u0 will become part of the code generated by FFC, and the complexity of u0 will be visible to FFC/UFL when estimating an integration order.
>> >
>>
>> With the above, the accuracy issue is shifted from interpolation to
>> integration. The functional above will be integrated approximately,
>> and one needs to make sure that the integration is sufficiently
>> accurate. If u0 is an Expression, the integration is exact but there
>> is an interpolation error.
>
> Another issue here is how the integral is evaluated. The important
> point is to not do the following:
>
> e = \int (u - u_h)^2 = \int u^2 + \int u_h^2 - \int 2*u*uh
>
> since this will lead to a subtraction of two numbers of approximately
> the same size (if the error is small).
>
> That is why errornorm interpolates u and u_h into the same space, then
> subtracts the dof values, and then performs the integration.

A better fix would be to make FFC not do that, since the
same situation can occur in other cases.

Martin

Revision history for this message
Anders Logg (logg) wrote :

On Mon, Aug 13, 2012 at 03:58:34PM -0000, Martin Sandve Alnæs wrote:
> On 13 August 2012 17:35, Anders Logg <email address hidden> wrote:
> > On Mon, Aug 13, 2012 at 12:27:23PM -0000, Garth Wells wrote:
> >> On 13 August 2012 12:35, Martin Sandve Alnæs <email address hidden> wrote:
> >> > In the email prior to this report, you say that:
> >> > """
> >> > 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)
> >> > """
> >> >
> >> > AFAIK this is _not_ correct. When making an Expression("..."), FFC does
> >> > not even know what the expression string is, and the degree that DOLFIN
> >> > assigns is independent of the expression string.
> >> >
> >>
> >> Yes, it is independent of the string but DOLFIN does not determine
> >> the degree. FFC makes a guess at the degree. The code is in
> >> ffc/ffc/analysis.py in the function
> >>
> >> _auto_select_degree(elements)
> >>
> >>
> >> > If, on the other hand, you use an UFL expression:
> >> > u0 = 1.0 + x[0]*sin(pi*x[0])
> >> > and compile a functional
> >> > error = (u-u0)**2*dx
> >> > then no interpolation of u0 takes place, since the exact expression for u0 will become part of the code generated by FFC, and the complexity of u0 will be visible to FFC/UFL when estimating an integration order.
> >> >
> >>
> >> With the above, the accuracy issue is shifted from interpolation to
> >> integration. The functional above will be integrated approximately,
> >> and one needs to make sure that the integration is sufficiently
> >> accurate. If u0 is an Expression, the integration is exact but there
> >> is an interpolation error.
> >
> > Another issue here is how the integral is evaluated. The important
> > point is to not do the following:
> >
> > e = \int (u - u_h)^2 = \int u^2 + \int u_h^2 - \int 2*u*uh
> >
> > since this will lead to a subtraction of two numbers of approximately
> > the same size (if the error is small).
> >
> > That is why errornorm interpolates u and u_h into the same space, then
> > subtracts the dof values, and then performs the integration.
>
> A better fix would be to make FFC not do that, since the
> same situation can occur in other cases.

Not do what? Expand sums?

The interpolation happens in DOLFIN and the errornorm function was
added just to avoid the type of numerical instabilities that may
result from subtracting large numbers.

--
Anders

Revision history for this message
Martin Sandve Alnæs (martinal) wrote :
Download full text (3.5 KiB)

On 14 August 2012 13:01, Anders Logg <email address hidden> wrote:
> On Mon, Aug 13, 2012 at 03:58:34PM -0000, Martin Sandve Alnæs wrote:
>> On 13 August 2012 17:35, Anders Logg <email address hidden> wrote:
>> > On Mon, Aug 13, 2012 at 12:27:23PM -0000, Garth Wells wrote:
>> >> On 13 August 2012 12:35, Martin Sandve Alnæs <email address hidden> wrote:
>> >> > In the email prior to this report, you say that:
>> >> > """
>> >> > 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)
>> >> > """
>> >> >
>> >> > AFAIK this is _not_ correct. When making an Expression("..."), FFC does
>> >> > not even know what the expression string is, and the degree that DOLFIN
>> >> > assigns is independent of the expression string.
>> >> >
>> >>
>> >> Yes, it is independent of the string but DOLFIN does not determine
>> >> the degree. FFC makes a guess at the degree. The code is in
>> >> ffc/ffc/analysis.py in the function
>> >>
>> >> _auto_select_degree(elements)
>> >>
>> >>
>> >> > If, on the other hand, you use an UFL expression:
>> >> > u0 = 1.0 + x[0]*sin(pi*x[0])
>> >> > and compile a functional
>> >> > error = (u-u0)**2*dx
>> >> > then no interpolation of u0 takes place, since the exact expression for u0 will become part of the code generated by FFC, and the complexity of u0 will be visible to FFC/UFL when estimating an integration order.
>> >> >
>> >>
>> >> With the above, the accuracy issue is shifted from interpolation to
>> >> integration. The functional above will be integrated approximately,
>> >> and one needs to make sure that the integration is sufficiently
>> >> accurate. If u0 is an Expression, the integration is exact but there
>> >> is an interpolation error.
>> >
>> > Another issue here is how the integral is evaluated. The important
>> > point is to not do the following:
>> >
>> > e = \int (u - u_h)^2 = \int u^2 + \int u_h^2 - \int 2*u*uh
>> >
>> > since this will lead to a subtraction of two numbers of approximately
>> > the same size (if the error is small).
>> >
>> > That is why errornorm interpolates u and u_h into the same space, then
>> > subtracts the dof values, and then performs the integration.
>>
>> A better fix would be to make FFC not do that, since the
>> same situation can occur in other cases.
>
> Not do what? Expand sums?

Expand "product of sums" into "sum of products".

> The interpolation happens in DOLFIN and the errornorm function was
> added just to avoid the type of numerical instabilities that may
> result from subtracting large numbers.

Yes, and the same nu...

Read more...

Revision history for this message
Garth Wells (garth-wells) wrote :
Download full text (3.8 KiB)

On 14 August 2012 12:01, Anders Logg <email address hidden> wrote:
> On Mon, Aug 13, 2012 at 03:58:34PM -0000, Martin Sandve Alnæs wrote:
>> On 13 August 2012 17:35, Anders Logg <email address hidden> wrote:
>> > On Mon, Aug 13, 2012 at 12:27:23PM -0000, Garth Wells wrote:
>> >> On 13 August 2012 12:35, Martin Sandve Alnæs <email address hidden> wrote:
>> >> > In the email prior to this report, you say that:
>> >> > """
>> >> > 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)
>> >> > """
>> >> >
>> >> > AFAIK this is _not_ correct. When making an Expression("..."), FFC does
>> >> > not even know what the expression string is, and the degree that DOLFIN
>> >> > assigns is independent of the expression string.
>> >> >
>> >>
>> >> Yes, it is independent of the string but DOLFIN does not determine
>> >> the degree. FFC makes a guess at the degree. The code is in
>> >> ffc/ffc/analysis.py in the function
>> >>
>> >> _auto_select_degree(elements)
>> >>
>> >>
>> >> > If, on the other hand, you use an UFL expression:
>> >> > u0 = 1.0 + x[0]*sin(pi*x[0])
>> >> > and compile a functional
>> >> > error = (u-u0)**2*dx
>> >> > then no interpolation of u0 takes place, since the exact expression for u0 will become part of the code generated by FFC, and the complexity of u0 will be visible to FFC/UFL when estimating an integration order.
>> >> >
>> >>
>> >> With the above, the accuracy issue is shifted from interpolation to
>> >> integration. The functional above will be integrated approximately,
>> >> and one needs to make sure that the integration is sufficiently
>> >> accurate. If u0 is an Expression, the integration is exact but there
>> >> is an interpolation error.
>> >
>> > Another issue here is how the integral is evaluated. The important
>> > point is to not do the following:
>> >
>> > e = \int (u - u_h)^2 = \int u^2 + \int u_h^2 - \int 2*u*uh
>> >
>> > since this will lead to a subtraction of two numbers of approximately
>> > the same size (if the error is small).
>> >
>> > That is why errornorm interpolates u and u_h into the same space, then
>> > subtracts the dof values, and then performs the integration.
>>
>> A better fix would be to make FFC not do that, since the
>> same situation can occur in other cases.
>
> Not do what? Expand sums?
>

I think Martin is referring to FFC automatically selecting polynomial
degrees for functions. I think that the present situation is the
lesser of two evils.

> The interpolation happens in DOLFIN and the errornorm function was
> added just to avoi...

Read more...

Revision history for this message
Martin Sandve Alnæs (martinal) wrote :
Download full text (3.4 KiB)

On 14 August 2012 13:31, Garth Wells <email address hidden> wrote:
> On 14 August 2012 12:01, Anders Logg <email address hidden> wrote:
>> On Mon, Aug 13, 2012 at 03:58:34PM -0000, Martin Sandve Alnæs wrote:
>>> On 13 August 2012 17:35, Anders Logg <email address hidden> wrote:
>>> > On Mon, Aug 13, 2012 at 12:27:23PM -0000, Garth Wells wrote:
>>> >> On 13 August 2012 12:35, Martin Sandve Alnæs <email address hidden> wrote:
>>> >> > In the email prior to this report, you say that:
>>> >> > """
>>> >> > 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)
>>> >> > """
>>> >> >
>>> >> > AFAIK this is _not_ correct. When making an Expression("..."), FFC does
>>> >> > not even know what the expression string is, and the degree that DOLFIN
>>> >> > assigns is independent of the expression string.
>>> >> >
>>> >>
>>> >> Yes, it is independent of the string but DOLFIN does not determine
>>> >> the degree. FFC makes a guess at the degree. The code is in
>>> >> ffc/ffc/analysis.py in the function
>>> >>
>>> >> _auto_select_degree(elements)
>>> >>
>>> >>
>>> >> > If, on the other hand, you use an UFL expression:
>>> >> > u0 = 1.0 + x[0]*sin(pi*x[0])
>>> >> > and compile a functional
>>> >> > error = (u-u0)**2*dx
>>> >> > then no interpolation of u0 takes place, since the exact expression for u0 will become part of the code generated by FFC, and the complexity of u0 will be visible to FFC/UFL when estimating an integration order.
>>> >> >
>>> >>
>>> >> With the above, the accuracy issue is shifted from interpolation to
>>> >> integration. The functional above will be integrated approximately,
>>> >> and one needs to make sure that the integration is sufficiently
>>> >> accurate. If u0 is an Expression, the integration is exact but there
>>> >> is an interpolation error.
>>> >
>>> > Another issue here is how the integral is evaluated. The important
>>> > point is to not do the following:
>>> >
>>> > e = \int (u - u_h)^2 = \int u^2 + \int u_h^2 - \int 2*u*uh
>>> >
>>> > since this will lead to a subtraction of two numbers of approximately
>>> > the same size (if the error is small).
>>> >
>>> > That is why errornorm interpolates u and u_h into the same space, then
>>> > subtracts the dof values, and then performs the integration.
>>>
>>> A better fix would be to make FFC not do that, since the
>>> same situation can occur in other cases.
>>
>> Not do what? Expand sums?
>>
>
> I think Martin is referring to FFC automatically selecting polynomial
> degrees for functions. I think th...

Read more...

Garth Wells (garth-wells)
Changed in dolfin:
status: New → In Progress
Garth Wells (garth-wells)
Changed in dolfin:
assignee: nobody → Garth Wells (garth-wells)
status: In Progress → Fix Committed
Marie Rognes (meg-simula)
Changed in dolfin:
status: Fix Committed → Fix Released
To post a comment you must log in.
This report contains Public information  
Everyone can see this information.
Subscribing...

Other bug subscribers

Related questions

Remote bug watches

Bug watches keep track of this bug in other bug trackers.