Functions/co-functions and vectors/co-vectors

[thomasgibson]

I'm raising an issue here to keep track of this ongoing problem we have with describing a function vs a co-function.

We know mathematically what the difference is: Functions are linear combinations of finite element basis functions. These are the usual firedrake.Functions (ufl.Coefficients). Cofunctions are bounded linear maps which take functions and returns some scalar.

I think this is what currently happens (please correct me if I'm wrong): Let f be some function. Then

g = assemble(inner(f, v)*dx)

for some test function v produces a Function, g, whose ith component is the evaulated inner product of f with the ith global finite element basis function. These are not functions! They are assembled functionals or co-functions are should be treated as such.

In Slate, there are AssembledVectors which are only used in local actions (mat-vec multiplication). In most cases, these are actual vectors formed from a finite element function, but there may be instances where these may be co-vectors corresponding to assembled 1-forms.

It makes sense to ask cofunctions for their arguments like a form, and ideally one should be able to use existing form operations on them with assembled or symbolic (UFL) forms. Basically, we need a way to distinguish functions/vectors and co-functions/co-vectors in Firedrake. This may require some changes to UFL.

So, questions:

What changes are required in UFL?
- New object: Co-coefficient?
- How should we handle expressions like this: let lv be a hypothetical co-function and f some coefficient. What do we do with:

lv + inner(f, v)*dx

Firedrake changes?
- Probably should add CoFunction or something like that.
- How to assemble expressions like the one above?
- Slate changes: AssembledVector/AssembledCoVector?

Please feel free to edit this post to add additional points/development discussions.

[danshapero]

I've run into this same problem in other contexts. An alternative to introducing CoFunction would be to make a class that represents the dual of a function space.

V = FunctionSpace(mesh, ...)
V_star = DualFunctionSpace(V)

f = Function(V)
v = TestFunction(V)
F = Function(V_star)
F.assign(assemble(f*v*dx))

x = Point(1.0, 1.0)
print(f(x))
>>>> 0.0

print(F(x))
>>>> Error! You can evaluate functions at point, not co-functions.

Not sure that it's any better or worse, just food for thought!

[danshapero]

I brought this up again in the slack channel the other day. Taking the dual of a function space might be a bad idea, because then you start having to worry about what happens when you take the dual twice. With that in mind, using a CoFunction class would be the smarter approach.

We also talked about how algebra on cofunctions would work -- what happens when you add lv to inner(f, v)*dx where lv is a CoFunction, f is a proper Function, and v is a test function? When you assemble a UFL bilinear form, you get a Matrix object. Trying to add a Matrix to another UFL form gives a type error:

mesh = ...
Q = FunctionSpace(mesh, ...)
u = TrialFunction(Q)
v = TestFunction(Q)
form = u * v * dx
M = assemble(form)
M + form
>>>> TypeError: unsupported operand type(s) for +: 'Matrix' and 'Form'

If you can't add a bilinear form and a matrix, then it's perfectly consistent if you also can't add a linear form and a cofunction.

[wence-]

Update, work in progress document describing the new algebraic constructs in UFL: https://www.overleaf.com/read/tqzqzdybsdsq

Comments to @dham.

[wence-]

Taking the dual of a function space might be a bad idea, because then you start having to worry about what happens when you take the dual twice.

Finite dimensional vector spaces are reflexive so V** is isometrically isomorphic to V, hence we can just identify V** as V.

[danshapero]

Taking the dual of a function space might be a bad idea, because then you start having to worry about what happens when you take the dual twice.

Finite dimensional vector spaces are reflexive so V** is isometrically isomorphic to V, hence we can just identify V** as V.

This came up in a slack discussion that isn't archived, but I think at the time someone objected on the grounds that L^1 isn't reflexive. I also had the impression from the discussions about Koki's work that, at the UFL level, a finite-dimensional basis hasn't been introduced yet. It seems like this is moving towards using CoFunctions rather than representing the dual spaces themselves, in which case the question about double-duals is a moot point. And either way I'm happy to see this moving forward!

[dham]

The plan is to have dual spaces. However, I don't think L^1 is a relevant space. Everything in UFL is in L^2 at least, and that space is reflexive (all Hilbert spaces are reflexive).

The objections to Koki's work was that it appeared to require UFL to know about the specific choice of basis. That UFL spaces are finite-dimensional and have a basis is, I believe, uncontroversial. The extensions I am proposing only assume the existence of a dual basis, not what it actually is. I think that won't cause too many fights.

[nbouziani]

I think this issue should be closed now that the dual space changes have been merged (see #2294)