Tidy up interpolation tsfc<->firedrake interface

[ReubenHill]

See issue #232. Coupled with firedrakeproject/firedrake#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.

tsfc/tsfc/driver.py

Line 136 in b043e86

builder.set_coordinates(mesh)

What does that function do?

tsfc/tsfc/kernel_interface/firedrake_loopy.py

Lines 221 to 229 in b043e86

	def set_coordinates(self, domain):
	"""Prepare the coordinate field.
	:arg domain: :class:`ufl.Domain`
	"""
	# Create a fake coordinate coefficient for a domain.
	f = Coefficient(FunctionSpace(domain, domain.ufl_coordinate_element()))
	self.domain_coordinate[domain] = f
	self.coordinates_arg = self._coefficient(f, "coords")

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?

[ReubenHill]

to_element should eventually become a UFL cofunction whenever they are implemented

[wence-]

TBH, this is a horrible variable name. I think what this means is that the kernel expects to be fed:

[mesh.coordinates] + extract_coefficients(expression)

Is that right?

[ReubenHill]

Didn't see this comment in time. Yes that is correct

[ReubenHill]

If you have a suggested better interface and/or variable name to say and I can make a new PR

[sv2518]

I think in various places in Firedrake we have similar variables to this e.g. for cell sizes. They are called something like needs_cell_sizes, maybe you just want to rename to needs_coords.

[ksagiyam]

The way coordinate is taken care of in the above is identical to how we do in KernelBuilder.set_coordinate, so maybe one potential future change is to make set_coordinate method also available in ExpressionKernelBuilder.

[ReubenHill]

I strongly considered doing this, but didn't want to make a "fake" version of that interface without recreating it more fully. Hopefully once your changes to that interface have landed this can be given another spin.

[ksagiyam]

Another thing I think we should think about next time is how to make ExpressionKernelBuilder more symmetric to KernelBuilder. I think instead of returning coefficients here, we should just call extract_coefficients(expr) again in Firedrake (this compares to form.coefficients() called in assemble.py) and split Functions if necessary when setting up ParLoop arguments (again, as is done in assemble.py).