Skip to content

Tidy up interpolation tsfc<->firedrake interface #233

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 3 commits into from
Nov 24, 2020
Merged

Tidy up interpolation tsfc<->firedrake interface #233

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
56c9b94
Remove coords argument in compile_expression_dual_evaluation
ReubenHill Oct 9, 2020
649a451
test compile_expression_dual_evaluation with symbolics
ReubenHill Oct 19, 2020
ed89971
Make clear if using manufactured coordinates coefficient
ReubenHill Nov 18, 2020
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
61 changes: 61 additions & 0 deletions tests/test_dual_evaluation.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,61 @@
import pytest
import ufl
from tsfc.finatinterface import create_element
from tsfc import compile_expression_dual_evaluation


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


def test_ufl_only_spatialcoordinate():
mesh = ufl.Mesh(ufl.VectorElement("P", ufl.triangle, 1))
V = ufl.FunctionSpace(mesh, ufl.FiniteElement("P", ufl.triangle, 2))
x, y = ufl.SpatialCoordinate(mesh)
expr = x*y - y**2 + x
W = V
to_element = create_element(W.ufl_element())
ast, oriented, needs_cell_sizes, coefficients, first_coeff_fake_coords, _ = compile_expression_dual_evaluation(expr, to_element, coffee=False)
assert first_coeff_fake_coords is True


def test_ufl_only_from_contravariant_piola():
mesh = ufl.Mesh(ufl.VectorElement("P", ufl.triangle, 1))
V = ufl.FunctionSpace(mesh, ufl.FiniteElement("RT", ufl.triangle, 1))
v = ufl.Coefficient(V)
expr = ufl.inner(v, v)
W = ufl.FunctionSpace(mesh, ufl.FiniteElement("P", ufl.triangle, 2))
to_element = create_element(W.ufl_element())
ast, oriented, needs_cell_sizes, coefficients, first_coeff_fake_coords, _ = compile_expression_dual_evaluation(expr, to_element, coffee=False)
assert first_coeff_fake_coords is True


def test_ufl_only_to_contravariant_piola():
mesh = ufl.Mesh(ufl.VectorElement("P", ufl.triangle, 1))
V = ufl.FunctionSpace(mesh, ufl.FiniteElement("P", ufl.triangle, 2))
v = ufl.Coefficient(V)
expr = ufl.as_vector([v, v])
W = ufl.FunctionSpace(mesh, ufl.FiniteElement("RT", ufl.triangle, 1))
to_element = create_element(W.ufl_element())
ast, oriented, needs_cell_sizes, coefficients, first_coeff_fake_coords, _ = compile_expression_dual_evaluation(expr, to_element, coffee=False)
assert first_coeff_fake_coords is True


def test_ufl_only_shape_mismatch():
mesh = ufl.Mesh(ufl.VectorElement("P", ufl.triangle, 1))
V = ufl.FunctionSpace(mesh, ufl.FiniteElement("RT", ufl.triangle, 1))
v = ufl.Coefficient(V)
expr = ufl.inner(v, v)
assert expr.ufl_shape == ()
W = V
to_element = create_element(W.ufl_element())
assert to_element.value_shape == (2,)
with pytest.raises(ValueError):
ast, oriented, needs_cell_sizes, coefficients, first_coeff_fake_coords, _ = compile_expression_dual_evaluation(expr, to_element, coffee=False)
27 changes: 15 additions & 12 deletions tsfc/driver.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -267,7 +267,7 @@ def name_multiindex(multiindex, name):
return builder.construct_kernel(kernel_name, impero_c, index_names, quad_rule)


def compile_expression_dual_evaluation(expression, to_element, coordinates, *,
def compile_expression_dual_evaluation(expression, to_element, *,
Copy link
Contributor Author

@ReubenHill ReubenHill Nov 24, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

to_element should eventually become a UFL cofunction whenever they are implemented

domain=None, interface=None,
parameters=None, coffee=False):
"""Compile a UFL expression to be evaluated against a compile-time known reference element's dual basis.
Expand All @@ -276,8 +276,7 @@ def compile_expression_dual_evaluation(expression, to_element, coordinates, *,

:arg expression: UFL expression
:arg to_element: A FInAT element for the target space
:arg coordinates: the coordinate function
:arg domain: optional UFL domain the expression is defined on (useful when expression contains no domain).
:arg domain: optional UFL domain the expression is defined on (required when expression contains no domain).
:arg interface: backend module for the kernel interface
:arg parameters: parameters object
:arg coffee: compile coffee kernel instead of loopy kernel
Expand Down Expand Up @@ -328,25 +327,29 @@ def compile_expression_dual_evaluation(expression, to_element, coordinates, *,
argument_multiindices = tuple(builder.create_element(arg.ufl_element()).get_indices()
for arg in arguments)

# Replace coordinates (if any)
domain = expression.ufl_domain()
if domain:
assert coordinates.ufl_domain() == domain
builder.domain_coordinate[domain] = coordinates
builder.set_cell_sizes(domain)
# Replace coordinates (if any) unless otherwise specified by kwarg
if domain is None:
domain = expression.ufl_domain()
assert domain is not None

# Collect required coefficients
first_coefficient_fake_coords = False
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
# Create a fake coordinate coefficient for a domain.
coords_coefficient = ufl.Coefficient(ufl.FunctionSpace(domain, domain.ufl_coordinate_element()))
builder.domain_coordinate[domain] = coords_coefficient
builder.set_cell_sizes(domain)
coefficients = [coords_coefficient] + coefficients
first_coefficient_fake_coords = True
builder.set_coefficients(coefficients)

Copy link
Contributor

@ksagiyam ksagiyam Nov 24, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Contributor Author

@ReubenHill ReubenHill Nov 24, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

# Split mixed coefficients
expression = ufl_utils.split_coefficients(expression, builder.coefficient_split)

# Translate to GEM
kernel_cfg = dict(interface=builder,
ufl_cell=coordinates.ufl_domain().ufl_cell(),
ufl_cell=domain.ufl_cell(),
argument_multiindices=argument_multiindices,
index_cache={},
scalar_type=parameters["scalar_type"])
Expand Down Expand Up @@ -431,7 +434,7 @@ def compile_expression_dual_evaluation(expression, to_element, coordinates, *,
# Handle kernel interface requirements
builder.register_requirements([ir])
# Build kernel tuple
return builder.construct_kernel(return_arg, impero_c, index_names)
return builder.construct_kernel(return_arg, impero_c, index_names, first_coefficient_fake_coords)


def lower_integral_type(fiat_cell, integral_type):
Expand Down
9 changes: 6 additions & 3 deletions tsfc/kernel_interface/firedrake_loopy.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,7 @@


# Expression kernel description type
ExpressionKernel = namedtuple('ExpressionKernel', ['ast', 'oriented', 'needs_cell_sizes', 'coefficients', 'tabulations'])
ExpressionKernel = namedtuple('ExpressionKernel', ['ast', 'oriented', 'needs_cell_sizes', 'coefficients', 'first_coefficient_fake_coords', 'tabulations'])


def make_builder(*args, **kwargs):
Expand Down Expand Up @@ -153,12 +153,14 @@ def register_requirements(self, ir):
provided by the kernel interface."""
self.oriented, self.cell_sizes, self.tabulations = check_requirements(ir)

def construct_kernel(self, return_arg, impero_c, index_names):
def construct_kernel(self, return_arg, impero_c, index_names, first_coefficient_fake_coords):
"""Constructs an :class:`ExpressionKernel`.

:arg return_arg: loopy.GlobalArg for the return value
:arg impero_c: gem.ImperoC object that represents the kernel
:arg index_names: pre-assigned index names
:arg first_coefficient_fake_coords: If true, the kernel's first
coefficient is a constructed UFL coordinate field
Copy link
Contributor

@wence- wence- Nov 24, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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?

Copy link
Contributor Author

@ReubenHill ReubenHill Nov 24, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Copy link
Contributor Author

@ReubenHill ReubenHill Nov 24, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Copy link
Contributor

@sv2518 sv2518 Nov 24, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

:returns: :class:`ExpressionKernel` object
"""
args = [return_arg]
Expand All @@ -173,7 +175,8 @@ def construct_kernel(self, return_arg, impero_c, index_names):
loopy_kernel = generate_loopy(impero_c, args, self.scalar_type,
"expression_kernel", index_names)
return ExpressionKernel(loopy_kernel, self.oriented, self.cell_sizes,
self.coefficients, self.tabulations)
self.coefficients, first_coefficient_fake_coords,
Copy link
Contributor

@ksagiyam ksagiyam Nov 24, 2020

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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).

self.tabulations)


class KernelBuilder(KernelBuilderBase):
Expand Down