Tidy up interpolation tsfc<->firedrake interface

[ReubenHill]

See issue #232. Coupled with https://github.com/firedrakeproject/firedrake/pull/1861. Please feel free to bikeshed this whilst it's WIP!

[ReubenHill]

So my broad impression is that we want to make ExpressionKernelBuilder more like KernelBuilder which seems to require much less firedrake specific data. Is that a fair interpretation? Any suggestions on how best to proceed @wence- ? I feel that what I've done so far kind of side steps the issue with inputting coordinates. They are still there just hidden within the domain.

[wence-]

So my broad impression is that we want to make ExpressionKernelBuilder more like KernelBuilder which seems to require much less firedrake specific data. Is that a fair interpretation? Any suggestions on how best to proceed @wence- ? I feel that what I've done so far kind of side steps the issue with inputting coordinates. They are still there just hidden within the domain.

Suppose I give you just the generated kernel (the loopy thing). And I say "this kernel matches up with interpolating this expression". And then the only thing I permit you to ask of the expression is "what are your coefficients" and "what is your domain". What information do I need to construct the parloop that executes the kernel?

[ReubenHill]

I would need to match the domain to a PyOP2 iteration set, match the coefficients to the par loop arguments (which we ought to do with something like coefficient_numbers that you get when you make a kernel from compiling a form), the access descriptors for the coefficients, and any special information about how the iterset should be traversed (I presume this is what oriented and needs_cell_sizes are for?)

Sound right or wrong?

Also this seems to be an answer to a different question - I am, as far as I'm aware, not in a situation where I have such a kernel - rather I'm trying to work out how you build one. The answer to your question seems one step ahead.

[wence-]

OK, I see. Basically, inside TSFC, if you want to disentangle the symbolics from the numerics you're only allowed to called UFL methods on the objects. IOW, you want to get to a situation where if you just had a fully symbolic expression, compile_expression_dual_evaluation would still work.

from ufl import *
from tsfc.finatinterface import create_element
V = VectorElement("P", triangle, 1)
mesh = Mesh(V)

V = FunctionSpace(mesh, FiniteElement("P", triangle, 2))
x, y = SpatialCoordinate(mesh)

expr = x*y - y**2 + x
element = create_element(V.ufl_element())
stuff = compile_expression_dual_evaluation(expr, element, coffee=False)

[ReubenHill]

OK, I see. Basically, inside TSFC, if you want to disentangle the symbolics from the numerics you're only allowed to called UFL methods on the objects. IOW, you want to get to a situation where if you just had a fully symbolic expression, compile_expression_dual_evaluation would still work.

from ufl import *
from tsfc.finatinterface import create_element
V = VectorElement("P", triangle, 1)
mesh = Mesh(V)

V = FunctionSpace(mesh, FiniteElement("P", triangle, 2))
x, y = SpatialCoordinate(mesh)

expr = x*y - y**2 + x
element = create_element(V.ufl_element())
stuff = compile_expression_dual_evaluation(expr, element, coffee=False)

This example should fail I believe - a function of spatial coordinates surely requires the coordinates data.

[wence-]

This example should fail I believe - a function of spatial coordinates surely requires the coordinates data.

Only when executing, not when compiling.

[ReubenHill]

This example should fail I believe - a function of spatial coordinates surely requires the coordinates data.

Only when executing, not when compiling.

So I can execute this

import ufl
from tsfc.finatinterface import create_element
from tsfc import compile_expression_dual_evaluation

def test_ufl_only_dependence():
    V = ufl.VectorElement("P", ufl.triangle, 1)
    mesh = ufl.Mesh(V)
    V = ufl.FunctionSpace(mesh, ufl.FiniteElement("P", ufl.triangle, 2))
    v = ufl.Coefficient(V)
    expr = ufl.inner(v, v)
    element = create_element(V.ufl_element())
    ast, oriented, needs_cell_sizes, coefficients, _ = compile_expression_dual_evaluation(expr, element, coffee=False)

and the ast has no code property. Am I aiming for similar behaviour with an expression of spatial coordinates?

[wence-]

and the ast has no code property. Am I aiming for similar behaviour with an expression of spatial coordinates?

I think we are talking at cross-purposes.

I claim that for generation of an interpolation kernel (which is what this is doing), I should need no global (Firedrake-level) information. But rather I should only need symbolic information (what the elements are, what the expression is, etc...).

I proposed one way of checking if this is the case right now is to feed purely symbolic objects (from UFL) into compile_expression_dual_evaluation. The kernel code that comes out should be identical to if I had fed in Firedrake objects.

[ReubenHill]

and the ast has no code property. Am I aiming for similar behaviour with an expression of spatial coordinates?

I think we are talking at cross-purposes.

I claim that for generation of an interpolation kernel (which is what this is doing), I should need no global (Firedrake-level) information. But rather I should only need symbolic information (what the elements are, what the expression is, etc...).

I proposed one way of checking if this is the case right now is to feed purely symbolic objects (from UFL) into compile_expression_dual_evaluation. The kernel code that comes out should be identical to if I had fed in Firedrake objects.

Okay - in a discussion with @dham he conjectured that if my expression is one of spatial coordinates or non-trivial pullback, I should need global domain information from the domain.coordinates coefficient (or we misunderstood each other). Why ought this not to be the case? Do the FInAT elements contain enough information already?

Put another way, how do I deal with has_type(expression, GeometricQuantity) == true to get past these lines:

    # Collect required coefficients
    coefficients = extract_coefficients(expression)
    if has_type(expression, GeometricQuantity) or any(fem.needs_coordinate_mapping(c.ufl_element()) for c in coefficients):
        coefficients = [coordinates] + coefficients
    builder.set_coefficients(coefficients)

(permalink)?

Does a correctly refactored ExpressionKernelBuilder not actually require the coordinates coefficient?

[wence-]

Does a correctly refactored ExpressionKernelBuilder not actually require the coordinates coefficient?

It needs a symbolic equivalent of that coefficient. A Coefficient(coordinates.ufl_function_space()) should work fine.

You shouldn't global domain information (numbers) when you're compiling a kernel. You only need the symbolic information that tells you how numbers are to be interpreted.

https://github.com/firedrakeproject/tsfc/blob/b043e8654a68345c524809953fc8f0705dc0d6eb/tsfc/driver.py#L136

What does that function do?

https://github.com/firedrakeproject/tsfc/blob/b043e8654a68345c524809953fc8f0705dc0d6eb/tsfc/kernel_interface/firedrake_loopy.py#L221-L229

So it makes a symbolic representation of the coordinate field from the mesh's coordinate element. No Firedrake object required.

[wence-]

Sorry, edited my comments for tetchiness.

[ReubenHill]

Okay so it looks like I've got everything compiling using UFL only. Now I need to work out how to get all the necessary coefficient processing information (which I need, do I need to split coefficients) out of compile_expression_dual_evaluation. I'm wary of doing much more until #234 has been reviewed since so much of what I will do will copy KernelBuilder.

[ReubenHill]

This now works - I wonder if there's a better way of dealing with the coordinates coefficient but perhaps this is good enough for now?