DOLFIN

inaccurate projection of Expressions onto high-order elements

Bug #1153105 reported by Andrew McRae
6
This bug affects 1 person
Affects Status Importance Assigned to Milestone
DOLFIN
Fix Released
Undecided
Unassigned
DOLFIN 1.2.0

Bug Description

*Edit, changed psiexpr to something more damning, updated errors below*

Code:
=========================================
from dolfin import *
mesh = UnitSquareMesh(16,16)
E = FunctionSpace(mesh, 'CG', 2)

psiexpr = Expression("x[0]*x[0] + 2*x[1]*x[1]")
psi1 = Function(interpolate(psiexpr,E))
psi2 = Function(project(psiexpr,E))

psi = TrialFunction(E)
chi = TestFunction(E)
a = psi*chi*dx
L = psiexpr*chi*dx
psi3 = Function(E)
solve(a==L,psi3)

print errornorm(psiexpr,psi1), "psi1" # interpolate
print errornorm(psiexpr,psi2), "psi2" # automatic project
print errornorm(psiexpr,psi3), "psi3" # manual project
print
print errornorm(psi1,psi2), "1 - 2 errornorm"
print errornorm(psi2,psi3), "2 - 3 errornorm"
print errornorm(psi3,psi1), "3 - 1 errornorm"
print
print sqrt(assemble((psi1-psi2)*(psi1-psi2)*dx)), "1 - 2 manual errornorm"
print sqrt(assemble((psi2-psi3)*(psi2-psi3)*dx)), "2 - 3 manual errornorm"
print sqrt(assemble((psi3-psi1)*(psi3-psi1)*dx)), "3 - 1 manual errornorm"
=========================================

Trunk output (dolfin version 1.1.0+, notice the exponents below):
1.27953563066e-15 psi1
0.000514738137008 psi2
0.000514742502223 psi3

0.000514738137008 1 - 2 errornorm
1.46884296622e-07 2 - 3 errornorm
0.000514742502223 3 - 1 errornorm

0.000514738137058 1 - 2 manual errornorm
1.46906436962e-07 2 - 3 manual errornorm
0.000514742502283 3 - 1 manual errornorm

FEniCS dev PPA output output (dolfin version 1.1.0):
1.27953563066e-15 psi1
3.6480070749e-07 psi2
1.35381390174e-15 psi3

3.64800707454e-07 1 - 2 errornorm
3.64800707458e-07 2 - 3 errornorm
3.44355978603e-16 3 - 1 errornorm

3.64805838291e-07 1 - 2 manual errornorm
3.64936234084e-07 2 - 3 manual errornorm
9.23118756652e-09 3 - 1 manual errornorm

FEniCS stable PPA output (dolfin version 1.1.0): same as Dev

Ubuntu-packaged FEniCS output (dolfin version 1.0.0):
1.2825848101e-15 psi1
3.39883963619e-07 psi2
1.37463683699e-15 psi3

3.39883963621e-07 1 - 2 errornorm
3.39883963643e-07 2 - 3 errornorm
3.85823947936e-16 3 - 1 errornorm

3.39911250342e-07 1 - 2 manual errornorm
3.40037905381e-07 2 - 3 manual errornorm
6.49275067829e-09 3 - 1 manual errornorm

Aside:
With dolfin 1.0.0, the program runs in well under a second. With all the others, it takes several seconds to run. Running cProfile shows the culprit:
   ncalls tottime percall cumtime percall filename:lineno(function)
        4 4.839 1.210 4.839 1.210 {_fem.new_DofMap}
        7 0.250 0.036 0.250 0.036 {_function.Function_interpolate}
<snip>

See original description
Revision history for this message
Andrew McRae (andymc) wrote :

FunctionSpaces DG-0, CG-1, DG-1 do not show this behaviour.

But CG-2, DG-2, CG-3, DG-3, etc. do.

summary: - possibly inaccurate FE solves in trunk
+ possibly inaccurate FE solves for high-order elements in trunk
Revision history for this message
Patrick Farrell (pefarrell) wrote : Re: possibly inaccurate FE solves for high-order elements in trunk

You could try doing a recursive bisection to see which revision caused the change. That is, if r1 is a known good revision (1.1.0 in your case) and r2 is a known bad revision, try the one in the middle, r3; if it fails, recursively apply the bisection to [r1, r3]; if it succeeds, recursively apply the bisection to [r2, r3].

Andrew McRae (andymc)
description: updated
Revision history for this message
Andrew McRae (andymc) wrote :

(Having lots of issues with putting together builds between 10th Jan and the start of March, so I'm giving up on the bisection)

Also, there seem to be two separate issues here:
1 [minor]) 'project' gives a noticeably worse result than manual projection. This is true for 'all' spaces in 1.1.0 and earlier, and spaces P1DG, P1, P0 in trunk.
Dumping fields to file makes this clear (e.g. the third entry should be exactly 1):

interpolate:
<dof index="0" value="0.9384765625" cell_index="30" cell_dof_index="5"/>
<dof index="1" value="1.001953125" cell_index="30" cell_dof_index="3"/>
<dof index="2" value="1" cell_index="30" cell_dof_index="1"/>
...
project:
<dof index="0" value="0.93847669325847582" cell_index="30" cell_dof_index="5"/>
<dof index="1" value="1.0019532509821611" cell_index="30" cell_dof_index="3"/>
<dof index="2" value="0.99999970541656547" cell_index="30" cell_dof_index="1"/>

manual project:
<dof index="0" value="0.93847656249999989" cell_index="30" cell_dof_index="5"/>
<dof index="1" value="1.0019531249999998" cell_index="30" cell_dof_index="3"/>
<dof index="2" value="1.0000000000000004" cell_index="30" cell_dof_index="1"/>

2 [MAJOR]) something bad happened to high-order function spaces (P2, P2DG, P3, ...) between 1.1.0 and trunk.

interpolate:
<dof index="0" value="0.9384765625" cell_index="30" cell_dof_index="5"/>
<dof index="1" value="1.001953125" cell_index="30" cell_dof_index="3"/>
<dof index="2" value="1" cell_index="30" cell_dof_index="1"/>

project:
<dof index="0" value="0.93822036670058628" cell_index="30" cell_dof_index="5"/>
<dof index="1" value="1.0019410650648681" cell_index="30" cell_dof_index="3"/>
<dof index="2" value="0.99857413162619046" cell_index="30" cell_dof_index="1"/>

manual project:
<dof index="0" value="0.93822023560294543" cell_index="30" cell_dof_index="5"/>
<dof index="1" value="1.0019409387279414" cell_index="30" cell_dof_index="3"/>
<dof index="2" value="0.99857442679001707" cell_index="30" cell_dof_index="1"/>

Also, it doesn't seem to be a problem with Expression("..."); I can replace psiexpr by psi1 when creating psi2 and psi3, and I get exactly the same result.

Revision history for this message
Garth Wells (garth-wells) wrote :

Try using

    solver_type="lu"

when you call project. The will eliminate any errors due to the Kyrlov solver.

Revision history for this message
Andrew McRae (andymc) wrote :

Indeed, that solves the minor issue!

Revision history for this message
Marie Rognes (meg-simula) wrote :

(1) Try specifying the degree for the Expression (help(Expression)):

psiexpr = Expression("x[0]*x[0] + 2*x[1]*x[1]", degree=2)

(2) Run instant-clean

(3) See if the results are more as expected

[(2) should not be necessary, but that might be a separate bug]

Revision history for this message
Marie Rognes (meg-simula) wrote :

See Bug #1153508 for a separate report for instant-clean issue.

Changed in dolfin:
status: New → Confirmed
Revision history for this message
Johan Hake (johan-hake) wrote : Re: [Bug 1153105] Re: possibly inaccurate FE solves for high-order elements in trunk

On 03/11/2013 09:53 AM, Marie Rognes wrote:
> (1) Try specifying the degree for the Expression (help(Expression)):
>
> psiexpr = Expression("x[0]*x[0] + 2*x[1]*x[1]", degree=2)
>
> (2) Run instant-clean
>
> (3) See if the results are more as expected
>
> [(2) should not be necessary, but that might be a separate bug]

On might get the impression that the generated C++ Expression relies on
the degree. That is not the case. Only the form which uses the
expression, should be recompiled. Which seems to be the case, btw:

from dolfin import *
mesh = UnitSquareMesh(2,2)
# Triggers JIT compilation
e = Expression("x[0]+x[1]", degree=1)

# Triggers JIT compilation
assemble(e*dx, mesh=mesh)

# No not trigger JIT compilation
e = Expression("x[0]+x[1]", degree=4)

# Triggers JIT compilation
assemble(e*dx, mesh=mesh)

Johan

Revision history for this message
Marie Rognes (meg-simula) wrote : Re: possibly inaccurate FE solves for high-order elements in trunk

@Johan: Agree, but the recompilation does not seem to be triggered in general cf. Bug #1153508

Revision history for this message
Johan Hake (johan-hake) wrote : Re: [Bug 1153105] Re: possibly inaccurate FE solves for high-order elements in trunk

On 03/11/2013 10:37 AM, Marie Rognes wrote:
> @Johan: Agree, but the recompilation does not seem to be triggered in
> general cf. Bug #1153508

Yeah, I saw that.

J

Revision history for this message
Andrew McRae (andymc) wrote : Re: possibly inaccurate FE solves for high-order elements in trunk

@Marie: your suggestion for specifying the degree seemed to work.

Also, my earlier comment is incorrect ("Also, it doesn't seem to be a problem with Expression("..."); I can replace psiexpr by psi1 when creating psi2 and psi3, and I get exactly the same result."). This only happened if ran the 'normal' version, then ran the modified version without instant-cleaning.

Andrew McRae (andymc)
summary: - possibly inaccurate FE solves for high-order elements in trunk
+ inaccurate projection of Expressions onto high-order elements
Revision history for this message
Andrew McRae (andymc) wrote :

This bug can be closed - it was a side-effect of [url=https://bugs.launchpad.net/ufl/+bug/1153508]Bug #1153508[/url], which has now been fixed.

Marie Rognes (meg-simula)
Changed in dolfin:
status: Confirmed → Fix Committed
Johannes Ring (johannr)
Changed in dolfin:
milestone: none → 1.2.0
Johannes Ring (johannr)
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

Remote bug watches

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