Dear Will, On 01/02/17 11:17, William Booker wrote:
Dear David and Lawrence,
I updated my DIV definition to
# Define discrete divergence operator def div_u(u, p): return (dot(u, grad(p)))*dx(domain=p.ufl_domain()) + (jump(p)*dot( conditional ( ge(dot(f,n)('+'), 0. ), u('+')*theta , u('-')*(1-theta)) , n('-')))*dS
Did I mistake how I was supposed to implement the direction conditional?
"/usr/not-backed-up/firedrake/2016-12-14/firedrake/src/tsfc/gem/gem.py", line 555, in __new__ assert set(multiindex) <= set(expression.free_indices) AssertionError
You found a bug in the form compiler. Fixed here: https://github.com/firedrakeproject/tsfc/pull/95 This will be merged soon, hopefully, and you can then firedrake-update to get the fix. I not in passing that your conditional expression doesn't do what you want. You want: u('+')*theta + u('-')*(1 - theta) to look the same under interchange of + and -, so you should write: dot(conditional(ge(dot(f, n)('+'), 0), u('+')*theta + u('-')*(1 - theta), u('-')*theta + u('+')*(1 - theta)), n('-')) That way, if the selection expression is true, you get the first guy, and if it is false (you're "looking from the other side") you get the second one. With the form compiler fix, and that change, you appear to conserve energy in parallel too. I note in passing that your code will be very slow because you call solve in a loop over-and-over. You should refactor it to use a LinearVariationalSolver. This notebook https://nbviewer.jupyter.org/github/firedrakeproject/firedrake/blob/master/d... describes how (and why) you need to do this. Lawrence