Subject: [Ffc] Quadrature optimisation bug

From: "Garth N. Wells" <[email-masked]>
Date: Mon, 01 Feb 2010 12:24:33 +0000

There seems to be a nasty bug in the quadrature optimisation. I can
compute a result, but it's wrong. I can try to look into it in more
detail later, but I'm using a combination of BDM and P0 functions which
might be enough of a hint for someone.

Garth

---

From: Marie Rognes <[email-masked]>
Date: Mon, 01 Feb 2010 22:20:14 +0100

> might be enough of a hint for someone.

Could you give a test script?

--
Marie

---

From: "Garth N. Wells" <[email-masked]>
Date: Mon, 01 Feb 2010 12:30:31 +0000

> Could you give a test script?

I'm working on it. The model is quite complicated so I'm trying to wind
it back.

Garth

---

From: Anders Logg <[email-masked]>
Date: Mon, 1 Feb 2010 13:32:58 +0100

> I'm working on it. The model is quite complicated so I'm trying to wind

Hurry up! I want to make a release today. Did you have to go looking
for bugs just now? ;-)

--
Anders

---

From: "Garth N. Wells" <[email-masked]>
Date: Mon, 01 Feb 2010 14:24:37 +0000

> Could you give a test script?

Below is a script which reproduces a bug when the quadrature
optimisation is turned on. The norm of vectors assembled with and
without using optimisation dffer.

Garth

from dolfin import *

mesh = UnitSquare(3, 3)
n = FacetNormal(mesh)

BDM = FunctionSpace(mesh, "Brezzi-Douglas-Marini", 1)
P0 = FunctionSpace(mesh, "Discontinuous Lagrange", 0)
mixed_space = MixedFunctionSpace([BDM, P0])

# Function spaces and functions
r  = TestFunction(P0)
U0 = Function(mixed_space)
u0, s0 = split(U0)

U0.vector()[:] = 1.0

un   = dot(u0, n)/2.0
F    = 1.0/(s0**2 + 1.0)
L = r*F*un*ds

param0 = {"representation": "quadrature", "optimize": False}
param1 = {"representation": "quadrature", "optimize": True}
print "Vec norm (no opt):   ", assemble(L,
                             form_compiler_parameters=param0).norm("l2")
print "Vec norm (with opt): ", assemble(L,
                             form_compiler_parameters=param1).norm("l2")

---

From: Marie Rognes <[email-masked]>
Date: Mon, 01 Feb 2010 15:35:27 +0100

> without using optimisation dffer.

Same problem arises with CG.

This looks like Harish's old bug where:

       un = dot(u0, n)/2.0

does not work, but

    un   = 0.5*dot(u0, n)

does.

(Just as a sidenote.)

--
Marie

---

From: "Garth N. Wells" <[email-masked]>
Date: Mon, 01 Feb 2010 14:40:05 +0000

> This looks like Harish's old bug where:

It's a strange one. I thought that Kristian had fixed it.

Interesting is that

  un   = dot(u0, n)/2.0
  F    = 1.0/(s0**2 + 1.0)

doesn't work, but

  un   = dot(u0, n)/2.0
  F    = s0**2

does work, as does

  un   = 0.5*dot(u0, n)
  F    = 1.0/(s0**2 + 1.0)

Garth

---

From: Kristian Oelgaard <[email-masked]>
Date: Mon, 1 Feb 2010 15:52:56 +0100 (CET)

> It's a strange one. I thought that Kristian had fixed it.

Strange, I would have guess that the  'F' expression in the first example was the problem, but that it works in the last example is a mystery. I'll have a look at it.

Kristian

> _______________________________________________
> Mailing list: https://launchpad.net/~ffc
> Post to     : [email-masked]
> Unsubscribe : https://launchpad.net/~ffc
> More help   : https://help.launchpad.net/ListHelp
>

---

From: "Garth N. Wells" <[email-masked]>
Date: Mon, 01 Feb 2010 14:53:34 +0000

Kristian Oelgaard wrote:
> Strange, I would have guess that the  'F' expression in the first
> example was the problem, but that it works in the last example is a
> mystery. I'll have a look at it.
> 

Thanks!

Garth