Projection of x is not accurate

[Martin Sandve Alnæs]

I've tested that projecting x works without the scaling bug that was just fixed, using dimensions 1,2,3 and both DG and CG from 0 to 3 degrees. I print the max and min values of the vector of the projection function, and the values are _close_ to 0 and 1 but not to machine precision. The script is below.

There's up to 2.7% error in the 3D case. Is the projection form integrated accurately enough? All but the DG0 function space should be capable of representing x exactly. Not sure if this is a dolfin or ffc bug.

from dolfin import *

def mcx(dim):
    if dim == 1:
        mesh = UnitInterval(20)
        cell = interval
        x = cell.x
    if dim == 2:
        mesh = UnitSquare(20, 20)
        cell = triangle
        x = cell.x[0]
    if dim == 3:
        mesh = UnitCube(20, 20, 20)
        cell = tetrahedron
        x = cell.x[0]
    return mesh, cell, x

for dim in range(1, 4):
    mesh, cell, x = mcx(dim)
    minval, maxval = 1.0, 0.0
    #print dim, "DG"
    for degree in range(3):
        V = FunctionSpace(mesh, "DG", degree)
        u = project(x, V)
        #print dim, degree, u.vector().min(), u.vector().max()
        minval = min(minval, u.vector().min())
        maxval = max(maxval, u.vector().max())
    #print dim, "CG"
    for degree in range(1, 4):
        V = FunctionSpace(mesh, "CG", degree)
        u = project(x, V)
        #print dim, degree, u.vector().min(), u.vector().max()
        minval = min(minval, u.vector().min())
        maxval = max(maxval, u.vector().max())
    print minval, maxval

[Garth Wells (garth-wells)]

There are two sources of error:

1) Inexact quadrature

2) The project function uses an iterative solver

[Martin Sandve Alnæs (martinal)]

Note that when replacing cell.x[0] with Expression("x[0]"), the accuracy of the projection is much better, around 1e-6, but not machine precision. I suspect the choice of quadrature rule is the source of this.

Here's a modified version of the above script doing this:

from dolfin import *

def mcx(dim):
    if dim == 1:
        mesh = UnitInterval(20)
        cell = interval
        x = cell.x
    if dim == 2:
        mesh = UnitSquare(20, 20)
        cell = triangle
        x = cell.x[0]
    if dim == 3:
        mesh = UnitCube(20, 20, 20)
        cell = tetrahedron
        x = cell.x[0]
    return mesh, cell, x

for dim in range(1, 4):
    mesh, cell, x = mcx(dim)
    x = Expression("x[0]")

    print dim, "DG"
    minval, maxval = 1.0, 0.0
    for degree in range(3):
        V = FunctionSpace(mesh, "DG", degree)
        u = project(x, V)
        print dim, degree, u.vector().min(), u.vector().max()
        minval = min(minval, u.vector().min())
        maxval = max(maxval, u.vector().max())
    print minval, maxval

    print dim, "CG"
    minval, maxval = 1.0, 0.0
    for degree in range(1, 4):
        V = FunctionSpace(mesh, "CG", degree)
        u = project(x, V)
        print dim, degree, u.vector().min(), u.vector().max()
        minval = min(minval, u.vector().min())
        maxval = max(maxval, u.vector().max())
    print minval, maxval

[Garth Wells (garth-wells)]

Using an LU solver and elevating the quadrature degree by one leads to machine zero error in all k > 0 test cases.

Using Constant reduced the error because FFC will use a P1 Lagrange function for x, which leads to exact quadrature.

[Martin Sandve Alnæs (martinal)]

The below script shows some strange behaviour with 'auto' quadrature degree and DG elements.

The script integrates du/dx over a unit square with u = the projection of 'pi x' into some space.

For DG spaces of degree 1 and 2 and several choices of quadrature orders, the result is either completely wrong or NaN.

Am I missing something here or should this be working better?
Also note that for CG elements the result of 'auto' quadrature degree is worse than 0.

from dolfin import *

# Projection which takes form compiler parameters
def myproject(expr, space, params):
    if 1: # Set to 0 to use the regular project and ignore params
        U = Function(space)
        u = TrialFunction(space)
        v = TestFunction(space)
        a = inner(u,v)*dx
        L = inner(expr,v)*dx
        problem = VariationalProblem(a, L, form_compiler_parameters=params)
        problem.solve(U)
    else:
        U = project(expr, space)
    return U

def test(fam, deg, qd):
    ffc_opt = { 'representation': 'quadrature', 'quadrature_degree': qd }
    n = 90
    mesh = UnitSquare(n,n)
    V = FunctionSpace(mesh, fam, deg)
    u = myproject(3.14*triangle.x[0], V, ffc_opt)
    a = u.dx(0)*dx
    dudx = assemble(u.dx(0)*dx, form_compiler_parameters=ffc_opt)
    print fam, deg, qd, dudx, norm(u.vector())

# Summary of the output is embedded in comments below:
print "Correct value is", 3.14
test("CG", 1, 0) # 3.14
test("CG", 1, 1) # 3.14
test("CG", 1, 2) # 3.14
test("CG", 1, 'auto') # 3.12 (worse than qd=0!)

test("DG", 2, 0) # NaN
test("DG", 2, 1) # NaN
test("DG", 2, 2) # NaN
test("DG", 2, 'auto') # 2.61 == 3.14/1.2

test("DG", 1, 0) # NaN
test("DG", 1, 1) # NaN
test("DG", 1, 2) # 3.14
test("DG", 1, 'auto') # almost zero

[Garth Wells (garth-wells)]

The test code should use an LU solver to separate out solver issues (and inaccuracies) from integration issues.

[Martin Sandve Alnæs (martinal)]

Using cpp.solve with "cg" and "default" instead of VariationalProblem.solve fixes the NaN problem. The numbers are still way off though:

Correct value is 3.14
CG 1 0 3.14000030019 165.429684591
CG 1 1 3.14000030019 165.429684591
CG 1 2 3.14000021752 165.429687084
CG 1 auto 3.12655064226 165.420523845
DG 2 0 2.34727023515e-10 290.790009255
DG 2 1 2.34727023515e-10 290.790009255
DG 2 2 3.92499865818 399.67114339
DG 2 auto 2.61666666667 565.213492026
DG 1 0 -3.8036831125e-08 399.652641011
DG 1 1 -3.8036831125e-08 399.652641011
DG 1 2 3.14 399.669087621
DG 1 auto -4.37743013134e-15 399.65264101

[Garth Wells (garth-wells)]

Use

  solver_type="lu"

when calling 'project'.

[Martin Sandve Alnæs (martinal)]

I've added form_compiler_parameters to project so I can set the quadrature degree (passed on to the two assemble calls).

This code, using "lu", still shows some problems:

def test2(qd):
    ffc_opt = { 'representation': 'quadrature', 'quadrature_degree': qd }
    mesh = UnitSquare(50, 50)
    V = FunctionSpace(mesh, "DG", 1)
    u = project(triangle.x[0], V, solver_type="lu", form_compiler_parameters=ffc_opt)
    a = u.dx(0)*dx
    M = assemble(u.dx(0)*dx, form_compiler_parameters=ffc_opt)
    uv = u.vector()[:]
    print uv.min(), uv.max(), M

test2(1)
test2(2)
test2(3)
test2('auto')

Output:

nan nan nan
-1.11173074327e-17 1.0 1.0
-9.52912065661e-18 1.0 1.0
0.00666666666667 0.993333333333 8.77744057992e-16

with qd = 2, 3, it looks good, with 1 or 'auto' we have nan or zero.

[Garth Wells (garth-wells)]

Some of these results are anomalies, and some are to be expected.

In all cases, if the quadrature order for the projection operation is increased, then the computed results are as expected.

For some cases, the prescribed quadrature scheme is not sufficient to integrate the projection, hence the NaN

[Anders Logg (logg)]

On Mon, May 30, 2011 at 05:37:17PM -0000, Garth Wells wrote:
> Some of these results are anomalies, and some are to be expected.
>
> In all cases, if the quadrature order for the projection operation is
> increased, then the computed results are as expected.
>
> For some cases, the prescribed quadrature scheme is not sufficient to
> integrate the projection, hence the NaN

Are any of the results caused by the FFC not detecting the right
quadrature rule to use when autodetecting the degree?

--
Anders

[Garth Wells (garth-wells)]

The NaN is not caused by FFC autodetecting the degree.

The numerical error in some cases is due to FFC not increasing the degree when the spatial coordinate is used as a coefficient function.

[Martin Sandve Alnæs (martinal)]

On 30 May 2011 21:02, Anders Logg <email address hidden> wrote:
> On Mon, May 30, 2011 at 05:37:17PM -0000, Garth Wells wrote:
>> Some of these results are anomalies, and some are to be expected.
>>
>> In all cases, if the quadrature order for the projection operation is
>> increased, then the computed results are as expected.
>>
>> For some cases, the prescribed quadrature scheme is not sufficient to
>> integrate the projection, hence the NaN
>
> Are any of the results caused by the FFC not detecting the right
> quadrature rule to use when autodetecting the degree?
>
> --
> Anders

If the autodetection of degree doesn't even cover projection I
would say it should be more conservative...

Martin

[Martin Sandve Alnæs (martinal)]

Ok, so NaN was because I manually specified too low degree.

So we're left with 'auto' which gives a big error:
0.00666666666667 0.993333333333 8.77744057992e-16

This is probably caused by ffc not considering
spatial coordinates in the expression. Fair enough.

So auto probably chooses 2 for the bilinear form
and 1 for the linear form, keeping the matrix nonsingular
but reducing the accuracy of the right hand side.

So the only bug here is that FFC could look at
spatial coordinates in the expression for better
autodetection, which would make project() of
expressions in PyDOLFIN more robust.
I'll close this bug and file a blueprint with this point.

Martin

On 30 May 2011 22:22, Garth Wells <email address hidden> wrote:
> The NaN is not caused by FFC autodetecting the degree.
>
> The numerical error in some cases is due to FFC not increasing the
> degree when the spatial coordinate is used as a coefficient function.
>
> --
> You received this bug notification because you are a member of DOLFIN
> Core Team, which is subscribed to DOLFIN.
> https://bugs.launchpad.net/bugs/785874
>
> Title:
>  Projection of x is not accurate
>
> Status in DOLFIN:
>  New
>
> Bug description:
>  I've tested that projecting x works without the scaling bug that was
>  just fixed, using dimensions 1,2,3 and both DG and CG from 0 to 3
>  degrees. I print the max and min values of the vector of the
>  projection function, and the values are _close_ to 0 and 1 but not to
>  machine precision. The script is below.
>
>  There's up to 2.7% error in the 3D case. Is the projection form
>  integrated accurately enough? All but the DG0 function space should be
>  capable of representing x exactly. Not sure if this is a dolfin or ffc
>  bug.
>
>
>  from dolfin import *
>
>  def mcx(dim):
>      if dim == 1:
>          mesh = UnitInterval(20)
>          cell = interval
>          x = cell.x
>      if dim == 2:
>          mesh = UnitSquare(20, 20)
>          cell = triangle
>          x = cell.x[0]
>      if dim == 3:
>          mesh = UnitCube(20, 20, 20)
>          cell = tetrahedron
>          x = cell.x[0]
>      return mesh, cell, x
>
>  for dim in range(1, 4):
>      mesh, cell, x = mcx(dim)
>      minval, maxval = 1.0, 0.0
>      #print dim, "DG"
>      for degree in range(3):
>          V = FunctionSpace(mesh, "DG", degree)
>          u = project(x, V)
>          #print dim, degree, u.vector().min(), u.vector().max()
>          minval = min(minval, u.vector().min())
>          maxval = max(maxval, u.vector().max())
>      #print dim, "CG"
>      for degree in range(1, 4):
>          V = FunctionSpace(mesh, "CG", degree)
>          u = project(x, V)
>          #print dim, degree, u.vector().min(), u.vector().max()
>          minval = min(minval, u.vector().min())
>          maxval = max(maxval, u.vector().max())
>      print minval, maxval
>

[Anders Logg (logg)]

Is there some UFL function that could be used instead to select the
degree (something that includes spatial coordinates)?

--
Anders

On Mon, May 30, 2011 at 08:44:15PM -0000, Martin Sandve Alnæs wrote:
> Ok, so NaN was because I manually specified too low degree.
>
> So we're left with 'auto' which gives a big error:
> 0.00666666666667 0.993333333333 8.77744057992e-16
>
> This is probably caused by ffc not considering
> spatial coordinates in the expression. Fair enough.
>
> So auto probably chooses 2 for the bilinear form
> and 1 for the linear form, keeping the matrix nonsingular
> but reducing the accuracy of the right hand side.
>
> So the only bug here is that FFC could look at
> spatial coordinates in the expression for better
> autodetection, which would make project() of
> expressions in PyDOLFIN more robust.
> I'll close this bug and file a blueprint with this point.
>
> Martin
>
> On 30 May 2011 22:22, Garth Wells <email address hidden> wrote:
> > The NaN is not caused by FFC autodetecting the degree.
> >
> > The numerical error in some cases is due to FFC not increasing the
> > degree when the spatial coordinate is used as a coefficient function.
> >
> >
> > Title:
> >  Projection of x is not accurate
> >
> > Status in DOLFIN:
> >  New
> >
> > Bug description:
> >  I've tested that projecting x works without the scaling bug that was
> >  just fixed, using dimensions 1,2,3 and both DG and CG from 0 to 3
> >  degrees. I print the max and min values of the vector of the
> >  projection function, and the values are _close_ to 0 and 1 but not to
> >  machine precision. The script is below.
> >
> >  There's up to 2.7% error in the 3D case. Is the projection form
> >  integrated accurately enough? All but the DG0 function space should be
> >  capable of representing x exactly. Not sure if this is a dolfin or ffc
> >  bug.
> >
> >
> >  from dolfin import *
> >
> >  def mcx(dim):
> >      if dim == 1:
> >          mesh = UnitInterval(20)
> >          cell = interval
> >          x = cell.x
> >      if dim == 2:
> >          mesh = UnitSquare(20, 20)
> >          cell = triangle
> >          x = cell.x[0]
> >      if dim == 3:
> >          mesh = UnitCube(20, 20, 20)
> >          cell = tetrahedron
> >          x = cell.x[0]
> >      return mesh, cell, x
> >
> >  for dim in range(1, 4):
> >      mesh, cell, x = mcx(dim)
> >      minval, maxval = 1.0, 0.0
> >      #print dim, "DG"
> >      for degree in range(3):
> >          V = FunctionSpace(mesh, "DG", degree)
> >          u = project(x, V)
> >          #print dim, degree, u.vector().min(), u.vector().max()
> >          minval = min(minval, u.vector().min())
> >          maxval = max(maxval, u.vector().max())
> >      #print dim, "CG"
> >      for degree in range(1, 4):
> >          V = FunctionSpace(mesh, "CG", degree)
> >          u = project(x, V)
> >          #print dim, degree, u.vector().min(), u.vector().max()
> >          minval = min(minval, u.vector().min())
> >          maxval = max(maxval, u.vector().max())
> >      print minval, maxval
> >
>

[Martin Sandve Alnæs (martinal)]

There's two, don't remember what they do:
  def estimate_max_polynomial_degree(e, default_degree=1):
  def estimate_total_polynomial_degree(e, default_degree=1):
in algorithms/transformations.py (should rather be in analysis.py I guess).

[Anders Logg (logg)]

On Mon, May 30, 2011 at 09:53:42PM -0000, Martin Sandve Alnæs wrote:
> There's two, don't remember what they do:
> def estimate_max_polynomial_degree(e, default_degree=1):
> def estimate_total_polynomial_degree(e, default_degree=1):
> in algorithms/transformations.py (should rather be in analysis.py I guess).
>
> ** Changed in: dolfin
> Status: New => Invalid

And those include spatial coordinates?

--
Anders

[Martin Sandve Alnæs (martinal)]

On 31 May 2011 00:24, Anders Logg <email address hidden> wrote:
> On Mon, May 30, 2011 at 09:53:42PM -0000, Martin Sandve Alnæs wrote:
>> There's two, don't remember what they do:
>>   def estimate_max_polynomial_degree(e, default_degree=1):
>>   def estimate_total_polynomial_degree(e, default_degree=1):
>> in algorithms/transformations.py (should rather be in analysis.py I guess).
>>
>> ** Changed in: dolfin
>>        Status: New => Invalid
>
> And those include spatial coordinates?

Turns out they didn't. Just checked the code.
But it was easy to add. I'm commiting changes
to estimate_total_polynomial_degree now which
incorporate the spatial degree. Maybe this should
be used for assembling rhs and functionals, while
looking at elements are enough for the bilinear form?

PyDOLFIN could do something like

d = estimate_total_polynomial_degree(expr)
d = max(d, 1)
d = min(d, 8)

to limit the degree to some reasonable range in cases such as
expr = sin(x**5)*cos(y**5)
which would lead to a degree of (5+2)+(5+2)=14 with the current heuristics.
Look at the code and tests in the last commit for more details, it's
quite short.

Martin

[Garth Wells (garth-wells)]

On 06/06/11 10:41, Martin Sandve Alnæs wrote:
> On 31 May 2011 00:24, Anders Logg <email address hidden> wrote:
>> On Mon, May 30, 2011 at 09:53:42PM -0000, Martin Sandve Alnæs wrote:
>>> There's two, don't remember what they do:
>>> def estimate_max_polynomial_degree(e, default_degree=1):
>>> def estimate_total_polynomial_degree(e, default_degree=1):
>>> in algorithms/transformations.py (should rather be in analysis.py I guess).
>>>
>>> ** Changed in: dolfin
>>> Status: New => Invalid
>>
>> And those include spatial coordinates?
>
> Turns out they didn't. Just checked the code.
> But it was easy to add. I'm commiting changes
> to estimate_total_polynomial_degree now which
> incorporate the spatial degree. Maybe this should
> be used for assembling rhs and functionals, while
> looking at elements are enough for the bilinear form?
>
> PyDOLFIN could do something like
>
> d = estimate_total_polynomial_degree(expr)
> d = max(d, 1)
> d = min(d, 8)
>
> to limit the degree to some reasonable range in cases such as
> expr = sin(x**5)*cos(y**5)
> which would lead to a degree of (5+2)+(5+2)=14 with the current heuristics.
> Look at the code and tests in the last commit for more details, it's
> quite short.
>

We have the same issue of order blow-out for problems with lots of
coefficients. I'm therefore inclined not to include any heuristics, and
leave it up to the user.

Garth

> Martin
>

[Anders Logg (logg)]

On Mon, Jun 06, 2011 at 09:54:12AM -0000, Garth Wells wrote:
> On 06/06/11 10:41, Martin Sandve Alnæs wrote:
> > On 31 May 2011 00:24, Anders Logg <email address hidden> wrote:
> >> On Mon, May 30, 2011 at 09:53:42PM -0000, Martin Sandve Alnæs wrote:
> >>> There's two, don't remember what they do:
> >>> def estimate_max_polynomial_degree(e, default_degree=1):
> >>> def estimate_total_polynomial_degree(e, default_degree=1):
> >>> in algorithms/transformations.py (should rather be in analysis.py I guess).
> >>>
> >>> ** Changed in: dolfin
> >>> Status: New => Invalid
> >>
> >> And those include spatial coordinates?
> >
> > Turns out they didn't. Just checked the code.
> > But it was easy to add. I'm commiting changes
> > to estimate_total_polynomial_degree now which
> > incorporate the spatial degree. Maybe this should
> > be used for assembling rhs and functionals, while
> > looking at elements are enough for the bilinear form?
> >
> > PyDOLFIN could do something like
> >
> > d = estimate_total_polynomial_degree(expr)
> > d = max(d, 1)
> > d = min(d, 8)
> >
> > to limit the degree to some reasonable range in cases such as
> > expr = sin(x**5)*cos(y**5)
> > which would lead to a degree of (5+2)+(5+2)=14 with the current heuristics.
> > Look at the code and tests in the last commit for more details, it's
> > quite short.
> >
>
> We have the same issue of order blow-out for problems with lots of
> coefficients. I'm therefore inclined not to include any heuristics, and
> leave it up to the user.

I think we should try to find some good heuristics. Setting "auto"
should try to do something sensible, but it needs to be well
documented what it does.

--
Anders

[Martin Sandve Alnæs (martinal)]

On 6 June 2011 11:54, Garth Wells <email address hidden> wrote:
> On 06/06/11 10:41, Martin Sandve Alnæs wrote:
>> On 31 May 2011 00:24, Anders Logg <email address hidden> wrote:
>>> On Mon, May 30, 2011 at 09:53:42PM -0000, Martin Sandve Alnæs wrote:
>>>> There's two, don't remember what they do:
>>>>   def estimate_max_polynomial_degree(e, default_degree=1):
>>>>   def estimate_total_polynomial_degree(e, default_degree=1):
>>>> in algorithms/transformations.py (should rather be in analysis.py I guess).
>>>>
>>>> ** Changed in: dolfin
>>>>        Status: New => Invalid
>>>
>>> And those include spatial coordinates?
>>
>> Turns out they didn't. Just checked the code.
>> But it was easy to add. I'm commiting changes
>> to estimate_total_polynomial_degree now which
>> incorporate the spatial degree. Maybe this should
>> be used for assembling rhs and functionals, while
>> looking at elements are enough for the bilinear form?
>>
>> PyDOLFIN could do something like
>>
>> d = estimate_total_polynomial_degree(expr)
>> d = max(d, 1)
>> d = min(d, 8)
>>
>> to limit the degree to some reasonable range in cases such as
>> expr = sin(x**5)*cos(y**5)
>> which would lead to a degree of (5+2)+(5+2)=14 with the current heuristics.
>> Look at the code and tests in the last commit for more details, it's
>> quite short.
>>
>
> We have the same issue of order blow-out for problems with lots of
> coefficients. I'm therefore inclined not to include any heuristics, and
> leave it up to the user.

Do you mean we should actually crash and burn with this line?
  f = assemble(sin(triangle.x[0]), mesh=mesh)
With the current heuristic this will give degree 3,
1 from x and +2 from sin.

Martin

[Garth Wells (garth-wells)]

On 06/06/11 11:05, Martin Sandve Alnæs wrote:
> On 6 June 2011 11:54, Garth Wells <email address hidden> wrote:
>> On 06/06/11 10:41, Martin Sandve Alnæs wrote:
>>> On 31 May 2011 00:24, Anders Logg <email address hidden> wrote:
>>>> On Mon, May 30, 2011 at 09:53:42PM -0000, Martin Sandve Alnæs wrote:
>>>>> There's two, don't remember what they do:
>>>>> def estimate_max_polynomial_degree(e, default_degree=1):
>>>>> def estimate_total_polynomial_degree(e, default_degree=1):
>>>>> in algorithms/transformations.py (should rather be in analysis.py I guess).
>>>>>
>>>>> ** Changed in: dolfin
>>>>> Status: New => Invalid
>>>>
>>>> And those include spatial coordinates?
>>>
>>> Turns out they didn't. Just checked the code.
>>> But it was easy to add. I'm commiting changes
>>> to estimate_total_polynomial_degree now which
>>> incorporate the spatial degree. Maybe this should
>>> be used for assembling rhs and functionals, while
>>> looking at elements are enough for the bilinear form?
>>>
>>> PyDOLFIN could do something like
>>>
>>> d = estimate_total_polynomial_degree(expr)
>>> d = max(d, 1)
>>> d = min(d, 8)
>>>
>>> to limit the degree to some reasonable range in cases such as
>>> expr = sin(x**5)*cos(y**5)
>>> which would lead to a degree of (5+2)+(5+2)=14 with the current heuristics.
>>> Look at the code and tests in the last commit for more details, it's
>>> quite short.
>>>
>>
>> We have the same issue of order blow-out for problems with lots of
>> coefficients. I'm therefore inclined not to include any heuristics, and
>> leave it up to the user.
>
> Do you mean we should actually crash and burn with this line?
> f = assemble(sin(triangle.x[0]), mesh=mesh)
> With the current heuristic this will give degree 3,
> 1 from x and +2 from sin.
>

We obviously need an approach for functions from non-polynomial spaces.
What I'm not inclined towards is arbitrary thresholds for integrating
polynomial products.

Garth

> Martin
>

[Anders Logg (logg)]

On Mon, Jun 06, 2011 at 10:17:27AM -0000, Garth Wells wrote:
> On 06/06/11 11:05, Martin Sandve Alnæs wrote:
> > On 6 June 2011 11:54, Garth Wells <email address hidden> wrote:
> >> On 06/06/11 10:41, Martin Sandve Alnæs wrote:
> >>> On 31 May 2011 00:24, Anders Logg <email address hidden> wrote:
> >>>> On Mon, May 30, 2011 at 09:53:42PM -0000, Martin Sandve Alnæs wrote:
> >>>>> There's two, don't remember what they do:
> >>>>> def estimate_max_polynomial_degree(e, default_degree=1):
> >>>>> def estimate_total_polynomial_degree(e, default_degree=1):
> >>>>> in algorithms/transformations.py (should rather be in analysis.py I guess).
> >>>>>
> >>>>> ** Changed in: dolfin
> >>>>> Status: New => Invalid
> >>>>
> >>>> And those include spatial coordinates?
> >>>
> >>> Turns out they didn't. Just checked the code.
> >>> But it was easy to add. I'm commiting changes
> >>> to estimate_total_polynomial_degree now which
> >>> incorporate the spatial degree. Maybe this should
> >>> be used for assembling rhs and functionals, while
> >>> looking at elements are enough for the bilinear form?
> >>>
> >>> PyDOLFIN could do something like
> >>>
> >>> d = estimate_total_polynomial_degree(expr)
> >>> d = max(d, 1)
> >>> d = min(d, 8)
> >>>
> >>> to limit the degree to some reasonable range in cases such as
> >>> expr = sin(x**5)*cos(y**5)
> >>> which would lead to a degree of (5+2)+(5+2)=14 with the current heuristics.
> >>> Look at the code and tests in the last commit for more details, it's
> >>> quite short.
> >>>
> >>
> >> We have the same issue of order blow-out for problems with lots of
> >> coefficients. I'm therefore inclined not to include any heuristics, and
> >> leave it up to the user.
> >
> > Do you mean we should actually crash and burn with this line?
> > f = assemble(sin(triangle.x[0]), mesh=mesh)
> > With the current heuristic this will give degree 3,
> > 1 from x and +2 from sin.
> >
>
> We obviously need an approach for functions from non-polynomial spaces.
> What I'm not inclined towards is arbitrary thresholds for integrating
> polynomial products.

What if we by default set that threshold to the maximal element degree
+ k, where k is say 1? It would be enough to retain expected
convergence.

--
Anders

[Martin Sandve Alnæs (martinal)]

On 6 June 2011 12:33, Anders Logg <email address hidden> wrote:
> On Mon, Jun 06, 2011 at 10:17:27AM -0000, Garth Wells wrote:
>> On 06/06/11 11:05, Martin Sandve Alnæs wrote:
>> > On 6 June 2011 11:54, Garth Wells <email address hidden> wrote:
>> >> On 06/06/11 10:41, Martin Sandve Alnæs wrote:
>> >>> On 31 May 2011 00:24, Anders Logg <email address hidden> wrote:
>> >>>> On Mon, May 30, 2011 at 09:53:42PM -0000, Martin Sandve Alnæs wrote:
>> >>>>> There's two, don't remember what they do:
>> >>>>>   def estimate_max_polynomial_degree(e, default_degree=1):
>> >>>>>   def estimate_total_polynomial_degree(e, default_degree=1):
>> >>>>> in algorithms/transformations.py (should rather be in analysis.py I guess).
>> >>>>>
>> >>>>> ** Changed in: dolfin
>> >>>>>        Status: New => Invalid
>> >>>>
>> >>>> And those include spatial coordinates?
>> >>>
>> >>> Turns out they didn't. Just checked the code.
>> >>> But it was easy to add. I'm commiting changes
>> >>> to estimate_total_polynomial_degree now which
>> >>> incorporate the spatial degree. Maybe this should
>> >>> be used for assembling rhs and functionals, while
>> >>> looking at elements are enough for the bilinear form?
>> >>>
>> >>> PyDOLFIN could do something like
>> >>>
>> >>> d = estimate_total_polynomial_degree(expr)
>> >>> d = max(d, 1)
>> >>> d = min(d, 8)
>> >>>
>> >>> to limit the degree to some reasonable range in cases such as
>> >>> expr = sin(x**5)*cos(y**5)
>> >>> which would lead to a degree of (5+2)+(5+2)=14 with the current heuristics.
>> >>> Look at the code and tests in the last commit for more details, it's
>> >>> quite short.
>> >>>
>> >>
>> >> We have the same issue of order blow-out for problems with lots of
>> >> coefficients. I'm therefore inclined not to include any heuristics, and
>> >> leave it up to the user.
>> >
>> > Do you mean we should actually crash and burn with this line?
>> >   f = assemble(sin(triangle.x[0]), mesh=mesh)
>> > With the current heuristic this will give degree 3,
>> > 1 from x and +2 from sin.
>> >
>>
>> We obviously need an approach for functions from non-polynomial spaces.
>> What I'm not inclined towards is arbitrary thresholds for integrating
>> polynomial products.
>
> What if we by default set that threshold to the maximal element degree
> + k, where k is say 1? It would be enough to retain expected
> convergence.
>
> --
> Anders

That's definitely not enough (e.g. L = x**8*v*dx) and would be basically
the same as throwing away the entire total polynomial degree algorithm.

It might be necessary/desirable to use different rules for bilinear
forms and linear forms/functionals?

Martin

[Garth Wells (garth-wells)]

This discussion should be moved to the mailing list or to a Blueprint.