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?

This gave me the following error:
  File "2d_compressible_stratified_2.py", line 180, in <module>
    solve(a == L, out, solver_parameters={'ksp_rtol': 1e-14})
  File "/usr/not-backed-up/firedrake/2016-12-14/firedrake/src/firedrake/firedrake/solving.py", line 122, in solve
    _solve_varproblem(*args, **kwargs)
  File "/usr/not-backed-up/firedrake/2016-12-14/firedrake/src/firedrake/firedrake/solving.py", line 151, in _solve_varproblem
    appctx=appctx)
  File "/usr/not-backed-up/firedrake/2016-12-14/firedrake/src/firedrake/firedrake/variational_solver.py", line 262, in __init__
    super(LinearVariationalSolver, self).__init__(*args, **kwargs)
  File "/usr/not-backed-up/firedrake/2016-12-14/firedrake/src/firedrake/firedrake/variational_solver.py", line 132, in __init__
    appctx=appctx)
  File "/usr/not-backed-up/firedrake/2016-12-14/firedrake/src/firedrake/firedrake/solving_utils.py", line 213, in __init__
    appctx=appctx)
  File "/usr/not-backed-up/firedrake/2016-12-14/firedrake/src/firedrake/firedrake/assemble.py", line 114, in allocate_matrix
    allocate_only=True)
  File "<decorator-gen-279>", line 2, in _assemble
  File "/usr/not-backed-up/firedrake/2016-12-14/firedrake/src/firedrake/firedrake/utils.py", line 62, in wrapper
    return f(*args, **kwargs)
  File "/usr/not-backed-up/firedrake/2016-12-14/firedrake/src/firedrake/firedrake/assemble.py", line 182, in _assemble
    inverse=inverse)
  File "/usr/not-backed-up/firedrake/2016-12-14/firedrake/src/firedrake/firedrake/tsfc_interface.py", line 189, in compile_form
    number_map).kernels
  File "/usr/not-backed-up/firedrake/2016-12-14/firedrake/src/PyOP2/pyop2/caching.py", line 200, in __new__
    obj = make_obj()
  File "/usr/not-backed-up/firedrake/2016-12-14/firedrake/src/PyOP2/pyop2/caching.py", line 190, in make_obj
    obj.__init__(*args, **kwargs)
  File "/usr/not-backed-up/firedrake/2016-12-14/firedrake/src/firedrake/firedrake/tsfc_interface.py", line 110, in __init__
    tree = tsfc_compile_form(form, prefix=name, parameters=parameters)
  File "/usr/not-backed-up/firedrake/2016-12-14/firedrake/src/tsfc/tsfc/driver.py", line 51, in compile_form
    kernels.append(compile_integral(integral_data, fd, prefix, parameters))
  File "/usr/not-backed-up/firedrake/2016-12-14/firedrake/src/tsfc/tsfc/driver.py", line 160, in compile_integral
    ir = fem.compile_ufl(integrand, interior_facet=interior_facet, **config)
  File "/usr/not-backed-up/firedrake/2016-12-14/firedrake/src/tsfc/tsfc/fem.py", line 380, in compile_ufl
    result = map_expr_dags(translator, expressions)
  File "/usr/not-backed-up/firedrake/2016-12-14/firedrake/src/ufl/ufl/corealg/map_dag.py", line 86, in map_expr_dags
    r = handlers[v._ufl_typecode_](v, *[vcache[u] for u in v.ufl_operands])
  File "/usr/not-backed-up/firedrake/2016-12-14/firedrake/src/tsfc/tsfc/ufl2gem.py", line 86, in conditional
    indices)
  File "/usr/not-backed-up/firedrake/2016-12-14/firedrake/src/tsfc/gem/gem.py", line 49, in __call__
    obj = super(NodeMeta, self).__call__(*args, **kwargs)
  File "/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

Thanks 
Will

From: firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> on behalf of Lawrence Mitchell <lawrence.mitchell@imperial.ac.uk>
Sent: 31 January 2017 17:44:40
To: firedrake@imperial.ac.uk
Subject: Re: [firedrake] Parallel issues and projection issues
 


On 31/01/17 12:49, William Booker wrote:
> Dear David,
>
> The choice of theta is meant to be arbitrary in the scheme, and we
> have had thoughts of looking into an optimal theta for implementations.
>
> So I would prefer to keep the asymmetry in the scheme.

So the point is that if theta = 0.5 is not optimal then it will be
*mesh dependent* (because if theta != 0.5 then the best choice is not
invariant under renumbering of the mesh).

> I'll implement the conditional form in my DIV definition, and get back
> to you about the results.
>
>
> Just a quick query, would the facet normal n not have to be assigned
> to - or + in that expression?

Yes, you want to take an (arbitrarily) restricted side of the expression:

so use:

dot(f, n)('+')


> Also, is there a simpler way to implement the projections that I have ?
> Going forward with a system that effectively doubles the amount of
> variables I have is a bit cumbersome.

So the physical variables are effectively diagnosed from the
variational derivatives?

But then the evolution equation says how the physical variable evolves
according to the derivatives?

Presumably then the elimination of one set, or the other, of these
variables requires that you invert some operator.  Is this operator
cell-local?

Cheers,

Lawrence