Hey guys I've been trying to implement a simple IP viscosity scheme for DG, but was getting weird results, which I think I've tracked down to a bug related to fd_bendy. If I run the following: from firedrake import * mesh2d = UnitSquareMesh(1,1) V = VectorFunctionSpace(mesh2d, "DG", 1) uv = Function(V) U_test = TestFunction(V) n = FacetNormal(mesh2d) F1 = dot(jump(U_test[0], n), jump(uv[0], n))*dS F2 = -dot(avg(grad(U_test[0])), jump(uv[0], n))*dS F=F1+F2 F0 = dot(jump(U_test[0],n)-avg(grad(U_test[0])), jump(uv[0], n))*dS M=assemble(derivative(F, uv)).M.values M1=assemble(derivative(F1, uv)).M.values M2=assemble(derivative(F2, uv)).M.values M0=assemble(derivative(F0, uv)).M.values I really expect M, M1+M2 and M0 to be all the same. On firedrake master (and in fact Dolfin) these are all the same indeed. But in my fd_bendy enviroment (fd_bendy for ffc and ufl, bendy_changes for firedrake) the values for M are different. Final piece of information is that if I use parameters['coffee']={}, everything is fine again. Other than that I'm stuck though, so I'd much appreciate any pointers into what I could check next Cheers Stephan
On 2 Oct 2015, at 13:10, Stephan Kramer <s.kramer@imperial.ac.uk> wrote:
Hey guys
I've been trying to implement a simple IP viscosity scheme for DG, but was getting weird results, which I think I've tracked down to a bug related to fd_bendy.
If I run the following:
from firedrake import *
mesh2d = UnitSquareMesh(1,1) V = VectorFunctionSpace(mesh2d, "DG", 1) uv = Function(V) U_test = TestFunction(V)
n = FacetNormal(mesh2d) F1 = dot(jump(U_test[0], n), jump(uv[0], n))*dS F2 = -dot(avg(grad(U_test[0])), jump(uv[0], n))*dS F=F1+F2 F0 = dot(jump(U_test[0],n)-avg(grad(U_test[0])), jump(uv[0], n))*dS
M=assemble(derivative(F, uv)).M.values M1=assemble(derivative(F1, uv)).M.values M2=assemble(derivative(F2, uv)).M.values M0=assemble(derivative(F0, uv)).M.values
I really expect M, M1+M2 and M0 to be all the same. On firedrake master (and in fact Dolfin) these are all the same indeed. But in my fd_bendy enviroment (fd_bendy for ffc and ufl, bendy_changes for firedrake) the values for M are different. Final piece of information is that if I use parameters['coffee']={}, everything is fine again. Other than that I'm stuck though, so I'd much appreciate any pointers into what I could check next
So the fact that things work correctly with coffee optimisations off is probably enough of a pointer. Can you try with the current proposed COFFEE pull request branch that Fabio has been working on, enhance-rewriter-mode3-logflops Locally, that appears to fix this problem. Lawrence
On 02/10/15 05:25, Lawrence Mitchell wrote:
On 2 Oct 2015, at 13:10, Stephan Kramer <s.kramer@imperial.ac.uk> wrote:
Hey guys
I've been trying to implement a simple IP viscosity scheme for DG, but was getting weird results, which I think I've tracked down to a bug related to fd_bendy.
If I run the following:
from firedrake import *
mesh2d = UnitSquareMesh(1,1) V = VectorFunctionSpace(mesh2d, "DG", 1) uv = Function(V) U_test = TestFunction(V)
n = FacetNormal(mesh2d) F1 = dot(jump(U_test[0], n), jump(uv[0], n))*dS F2 = -dot(avg(grad(U_test[0])), jump(uv[0], n))*dS F=F1+F2 F0 = dot(jump(U_test[0],n)-avg(grad(U_test[0])), jump(uv[0], n))*dS
M=assemble(derivative(F, uv)).M.values M1=assemble(derivative(F1, uv)).M.values M2=assemble(derivative(F2, uv)).M.values M0=assemble(derivative(F0, uv)).M.values
I really expect M, M1+M2 and M0 to be all the same. On firedrake master (and in fact Dolfin) these are all the same indeed. But in my fd_bendy enviroment (fd_bendy for ffc and ufl, bendy_changes for firedrake) the values for M are different. Final piece of information is that if I use parameters['coffee']={}, everything is fine again. Other than that I'm stuck though, so I'd much appreciate any pointers into what I could check next
So the fact that things work correctly with coffee optimisations off is probably enough of a pointer. Can you try with the current proposed COFFEE pull request branch that Fabio has been working on, enhance-rewriter-mode3-logflops
What's that magic! That fixes things for me indeed! Any idea what this could have caused? Or is this just a case of code refactoring accidentally fixing a bug? Cheers Stephan
Locally, that appears to fix this problem.
Lawrence
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
On 2 Oct 2015, at 13:40, Stephan Kramer <s.kramer@imperial.ac.uk> wrote:
On 02/10/15 05:25, Lawrence Mitchell wrote:
On 2 Oct 2015, at 13:10, Stephan Kramer <s.kramer@imperial.ac.uk> wrote:
Hey guys
I've been trying to implement a simple IP viscosity scheme for DG, but was getting weird results, which I think I've tracked down to a bug related to fd_bendy.
If I run the following:
from firedrake import *
mesh2d = UnitSquareMesh(1,1) V = VectorFunctionSpace(mesh2d, "DG", 1) uv = Function(V) U_test = TestFunction(V)
n = FacetNormal(mesh2d) F1 = dot(jump(U_test[0], n), jump(uv[0], n))*dS F2 = -dot(avg(grad(U_test[0])), jump(uv[0], n))*dS F=F1+F2 F0 = dot(jump(U_test[0],n)-avg(grad(U_test[0])), jump(uv[0], n))*dS
M=assemble(derivative(F, uv)).M.values M1=assemble(derivative(F1, uv)).M.values M2=assemble(derivative(F2, uv)).M.values M0=assemble(derivative(F0, uv)).M.values
I really expect M, M1+M2 and M0 to be all the same. On firedrake master (and in fact Dolfin) these are all the same indeed. But in my fd_bendy enviroment (fd_bendy for ffc and ufl, bendy_changes for firedrake) the values for M are different. Final piece of information is that if I use parameters['coffee']={}, everything is fine again. Other than that I'm stuck though, so I'd much appreciate any pointers into what I could check next
So the fact that things work correctly with coffee optimisations off is probably enough of a pointer. Can you try with the current proposed COFFEE pull request branch that Fabio has been working on, enhance-rewriter-mode3-logflops
What's that magic! That fixes things for me indeed! Any idea what this could have caused? Or is this just a case of code refactoring accidentally fixing a bug?
This is a pretty large scale rewrite of coffee's rewrite engine. Adding new features, but also possibly making things more robust. You're welcome to try and find the culprit! But note there are about 75 commits here: https://github.com/coneoproject/COFFEE/pull/61 Lawrence
On 2 Oct 2015, at 13:10, Stephan Kramer <s.kramer@imperial.ac.uk> wrote:
Hey guys
I've been trying to implement a simple IP viscosity scheme for DG, but was getting weird results, which I think I've tracked down to a bug related to fd_bendy.
If I run the following:
from firedrake import *
mesh2d = UnitSquareMesh(1,1) V = VectorFunctionSpace(mesh2d, "DG", 1) uv = Function(V) U_test = TestFunction(V)
n = FacetNormal(mesh2d) F1 = dot(jump(U_test[0], n), jump(uv[0], n))*dS F2 = -dot(avg(grad(U_test[0])), jump(uv[0], n))*dS F=F1+F2 F0 = dot(jump(U_test[0],n)-avg(grad(U_test[0])), jump(uv[0], n))*dS
M=assemble(derivative(F, uv)).M.values M1=assemble(derivative(F1, uv)).M.values M2=assemble(derivative(F2, uv)).M.values M0=assemble(derivative(F0, uv)).M.values
Do you mind wrapping this up as a test (to merge against current master) that ensures we don't regress when merging bendy? Lawrence
More weirdness, I'm afraid. If I change the function space V to: V = FunctionSpace(mesh2d, "RT", 1) I still get different answers between M and M1+M2. This is using the bendy branches for firedrake/ffc/ufl and enhance-rewriter-mode3-logflops for COFFEE. Again M==M1+M2 using parameters['coffee']={} and I get the same answer for M in "flat" firedrake. Any other magic branches I could try? ;-) Seriously though, any suggestions are welcome Cheers Stephan On 02/10/15 13:10, Stephan Kramer wrote:
Hey guys
I've been trying to implement a simple IP viscosity scheme for DG, but was getting weird results, which I think I've tracked down to a bug related to fd_bendy.
If I run the following:
from firedrake import *
mesh2d = UnitSquareMesh(1,1) V = VectorFunctionSpace(mesh2d, "DG", 1) uv = Function(V) U_test = TestFunction(V)
n = FacetNormal(mesh2d) F1 = dot(jump(U_test[0], n), jump(uv[0], n))*dS F2 = -dot(avg(grad(U_test[0])), jump(uv[0], n))*dS F=F1+F2 F0 = dot(jump(U_test[0],n)-avg(grad(U_test[0])), jump(uv[0], n))*dS
M=assemble(derivative(F, uv)).M.values M1=assemble(derivative(F1, uv)).M.values M2=assemble(derivative(F2, uv)).M.values M0=assemble(derivative(F0, uv)).M.values
I really expect M, M1+M2 and M0 to be all the same. On firedrake master (and in fact Dolfin) these are all the same indeed. But in my fd_bendy enviroment (fd_bendy for ffc and ufl, bendy_changes for firedrake) the values for M are different. Final piece of information is that if I use parameters['coffee']={}, everything is fine again. Other than that I'm stuck though, so I'd much appreciate any pointers into what I could check next
Cheers Stephan
_______________________________________________ firedrake mailing list firedrake@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/firedrake
Hey guys More weirdness related to facet integrals from IP viscosity that only appears with COFFEE optimisations. This time both on the bendy branches and on master (of last week). In the code below, the matrices M, Ms and M1+M2 should all be the same. For DG1 they are not equal on a recent (last week, before any sprint merges) install of firedrake/ffc/ufl/coffee on master. It is also incorrect on an older install I had, so it's not a recent thing. DG1 does however produce the right answer on the bendy branches: firedrake (bendy_changes) + ffc/ufl (fd_bendy) + coffee (master). For RT1 on the other hand it's the opposite: the matrices are correct and equal on master, but with bendy some of them are different from the correct answer. All of this only if use COFFEE optimisation, i.e. with parameters['coffee']={} the matrices are correct and all equal. I realize at this point the best thing to do is probably to wait if this all gets magically fixed when the sprint is over - will report back if not - but just thought to let you guys know already Cheers Stephan from firedrake import * mesh2d = UnitSquareMesh(1,1) U = VectorFunctionSpace(mesh2d, "DG", 1) v = TestFunction(U) u = Function(U) #parameters['coffee']={} n = FacetNormal(mesh2d) def outer_jump(v, n): return outer(v('+'), n('+'))+outer(v('-'), n('-')) F1 = inner(outer_jump(v, n), outer_jump(u,n))*dS F2 = inner(outer_jump(v, n), outer_jump(n,u))*dS F = inner(outer_jump(v, n), outer_jump(u,n)+outer_jump(n,u))*dS M1 = assemble(derivative(F1, u)).M.values M2 = assemble(derivative(F2, u)).M.values Ms = assemble(derivative(F1+F2, u)).M.values M = assemble(derivative(F, u)).M.values print abs(M-(M1+M2)).max() print abs(M-Ms).max()
participants (2)
-
Lawrence Mitchell
-
Stephan Kramer