aa r X i v : . [ c s . M S ] J a n UFL Dual Spaces, a proposal
David HamJanuary 2021
Abstract
This white paper highlights current limitations in the algebraic closure Unified FormLanguage (UFL). UFL currently represents forms over finite element spaces, however finiteelement problems naturally result in objects in the dual to a finite element space, and op-erators mapping between primal and dual finite element spaces. This document sketchesthe relevant mathematical areas and proposes changes to the UFL language to support dualspaces as first class types in UFL.
UFL is the Unified Form Language. It is a language used to describe multilinear forms. Typi-cally these forms are to be used in a finite-dimensional variational problem, specifically a finiteelement problem. UFL has not, hitherto, concerned itself with the spaces to which forms be-long or into which they map: this has been the concern of the assembler, such as Firedrake orFEniCS and, to some extent, the linear algebra backend (PETSc in the relevant cases, thoughother choices would be possible). These spaces involve the duals to the spaces in which UFLobjects such as Coefficients and Arguments are defined. The absence of dual spaces in UFLmeans that the language is not closed: there are operations visible to the UFL user which resultin objects that are not UFL objects, or which rely on bespoke functions where a more closedsystem would have a first class operator. Examples of this include:1. Interpolation, which is a form whose range space is the dual of the space being interpo-lated into.2. Nonlinear constitutive laws for materials such as ice. In these cases, UFL-based systemscurrently rely on an unnecessary Galerkin projection to represent an operation which isnot a differential operator.3. The adjoint, which currently has to rely on intricate reasoning to work out which objectsare primal or dual, where this should be obvious from their type.This document is an outline proposal for the inclusion of dual spaces and related objects inUFL. It has particular implications for ongoing work in Firedrake on interpolation and opera-tors external to the UFL system.
We will assume that there is a domain Ω ⊂ R n for some n ≥ Definition 2.1 (Function space) A function space is a vector space over a field K comprising func-tions Ω → K m . V is a function space and f ∈ V then f associates a value in K m with each pointin the domain Ω . The field K will be either R or C . Definition 2.2 (linear k -form) A linear k -form is a linear mapping from k ≥ vector spaces to K.That is, a linear k-form is a function:V k − × V k − × . . . × V → K (1) for suitable vector spaces V . . . V k − . Note 2.1
The argument numbering in definition 2.2 is deliberately from right to left. This reflects thenumbering of arguments in UFL. It also caters for a difference in notation between forms and linearalgebra: the assembled tensors are the transpose of the form, so the function spaces number from left toright when we consider the eventual linear algebra operations.
Definition 2.3 (sesquilinear k -form) A sesquilinear k -form is a mapping from k ≥ vector spacesto K which is conjugate linear (also termed antilinear) in the final argument (for k > ), and linear inthe others. Definition 2.4 (Inner product)
An inner product defined on a vector space V is a sesquilinear 2-formV × V → K. Which is Hermitian symmetric and positive definite. That is, for any u , v ∈ V1. h u , v i = h v , u i h u , u i ≥ Definition 2.5 (Hilbert space)
A Hilbert space is a complete normed vector space equipped with aninner product. All of the function spaces we are interested in are Hilbert spaces.
Definition 2.6 ( L inner product) The L inner product between two functions u , v : Ω → K isgiven by: h u , v i L = Z Ω u · v d x (2) Definition 2.7 ( L ) The function space L ( Ω ) n : Ω → K n is the Hilbert space of functions where thenorm induced by the L inner product is finite:u ∈ L ( Ω ) n ⇔ q h u , u i L ∈ R (3) Definition 2.8 (Dual space)
If V is a vector space then V ∗ , the dual space to V, is the space of boundedlinear functionals V → K. For the finite-dimensional spaces with which we are concerned, it is easy to establish the fol-lowing:
Theorem 2.1
If V is a finite-dimensional vector space, then | V | = | V ∗ | Theorem 2.2 (Dual basis) If { φ i } is a basis for a finite vector space V then there exists { φ ∗ i } a basisfor V ∗ such that: φ ∗ i ( φ j ) = δ ij (4) The basis { φ ∗ i } is termed the dual basis . Where it is necessary to make the distinction, we will refer to the space to which a dual spaceis dual as the primal space and its basis as the primal basis .2 efinition 2.9 Reflexive vector space A reflexive vector space is one for which ( V ∗ ) ∗ is isomorphic toV under the canonical map. That is, we can identify ( V ∗ ) ∗ and V. Theorem 2.3
All finite-dimensional vector spaces are reflexive.
This is a completely intuitive result. If u ∈ V and v ∈ V ∗ then this induces a u ′ in ( V ∗ ) ∗ givenby: u ′ ( v ) = v ( u ) (5)We simply identify u ′ and u . Since the dimensions of the spaces match and the functionals arelinear, this identification is an isometric isomorphism. In other words we can’t distinguish thespaces so the identification is valid. Theorem 2.4
All Hilbert spaces are reflexive.
In definition 2.2 and elsewhere, we use a notation for multilinear forms which takes a Cartesianproduct of spaces and returns a scalar. For example for a 2-form: V × V → K (6)This notation emphasises that the form is a function which takes two vectors and returns ascalar. However we could just as well Curry this form and write: V → V → K (7)This notation is left-associative, so this means: V → ( V → K ) (8)In other words, a 2-form is a function which takes a vector and returns a 1-form. k -forms In this section we restrict our attention to real-valued function spaces so as not to have toconsider where all the conjugates should be placed. We will come back to the complex valuedcase.It follows immediately from definitions 2.2 and 2.8 that the space of linear 1-forms on vectorspace V is the same thing as V ∗ . This is nothing else than writing: V ∗ = V → R (9)A direct consequence of this is the following equivalence: V × V → R = V → V ∗ (10)This says that every form can equivalently be thought of as an operator mapping into the dualof its 0-number space. 3 The interpolation operator
Usually we are concerned with vector spaces which are finite-dimensional subspaces of larger,typically infinite dimensional, spaces. In the finite element method, we often define the basis ofthese finite spaces by defining the space and choosing the dual basis. The primal basis can thenbe derived by solving a linear system (see, for example, the finite element course). Let’s brieflyexamine the dual basis functions, also known as nodes. These are functionals that take in afunction and (linearly) produce a scalar. What do such objects look like? The simplest case isprobably the point evaluation functional, or Dirac delta distribution. For a function v : Ω → R ,the point evaluation functional at the point x is given by: δ x ( v ) = v ( x ) (11)Other common functionals include the integral of a function over one mesh cell, or the com-ponent of a vector value in a particular direction at a particular point. All of these have incommon that they are trivially extensible to a much larger space of functions than the finite ele-ment function space whose dual they define. For example, point evaluation functionals can beevaluated for any space of functions whose members are well-defined at the points in question.This is useful because int is a very common requirement to approximate some function froma more general space in a finite element function. This most commonly occurs when a forcingfunction or reference solution is known either from external (physical) data, or from an analyticexpression. Suppose v ∈ V and that { φ i } is a basis for V , with { φ ∗ i } the corresponding dualbasis. The for a function f : Ω → R such that the dual functionals are defined, the interpolationof f into v is given by: ˆ f = ˆ f i φ i (12)where: ˆ f i = φ ∗ i ( f ) (13)In other words, the i -th basis coefficient is the evaluation of the i -th dual basis functional on thefunction to be interpolated. We can call this operator I V ( f ) . Note 5.1
If f ∈ V then it is straightforward to show that I V ( f ) = f . That is to say, interpolation is aprojection. In fact, the interpolation operator is nothing other than the evaluation of the dual basis. Wecould instead define an operator, a generalisation of the delta distribution, which takes a func-tional v ∈ V ∗ and applies it to a function u ∈ U where U is a function space on which thefunctionals in V ∗ are defined. That is to say: δ ( u , v ) = v ( u ) (14)Notice that δ : U × V ∗ → R . Since V ∗ is itself a vector space, this makes δ a 2-form, albeit aslightly unusual one. We would be more used to forms in U × V → R . UFL is the Unified Form Language (Alnæs et al. 2014) also described in Homolya et al. 2018. Asone might expect, it provides support for k -linear forms. UFL is a symbolic algebra language.In common with other such languages, it provides for symbols to stand for unknown values,including the arguments to k -linear forms. UFL is a purely symbolic language which problemsolving environments such as Firedrake and FEniCS extend by subclassing in order to attachnumerical data and actually solve problems. In understanding how UFL objects do (or should)4ehave, it is sometimes necessary to consider how they would be integrated into the problemsolving environment, and especially what should happen when a form is assembled.In an attempt to avoid replicating the documentation of UFL, we will start here from UFL FunctionSpace and gloss over the cells, meshes, and finite elements from which these are con-structed.If V is a ufl.FunctionSpace then c = ufl.Coefficient(V) defines a known function in V .Conversely, a = ufl.Argument(V, 0) defines a placeholder symbol for an unknown functionin V . Consequently, f_0 = c * dx is a symbolic expression for the integral of c over the domain and is a real value. UFL callsthis a 0-form. In contrast, f_1 = a * dx represents the integration of the unknown function a over the domain. It’s therefore a 1-form, or a function in V → R .In order to make the argument order clear, we introduce a second function space, W . Wedefine b = ufl.Argument(W, 1) and create the simplest 2-form: f_2 = c * b * dx This is a function in W × V → R . Notice that, by convention, argument 1 comes beforeargument 0. The problem solving environment takes symbolic UFL forms and evaluates them. In Firedrakeand FEniCS, the operation which does this is called assembly. The problem solving environ-ment encodes in its subclass of ufl.FunctionSpace a basis for the function space. We willwrite φ i for the basis for V and ψ i for the basis for W . The problem solving environment willalso implement its own subclass of ufl.Coefficient , in Firedrake and FEniCS this is called a Function . If c is a Coefficient , then mathematically: c = c i φ i (15)The consequence of this is that the core of the implementation of a Function is a reference to a
FunctionSpace (which encodes { φ i } ) and a vector of values c i . The latter give the Function itsnumerical value, while the former gives that value a mathematical meaning.If we execute: g = assemble(f_0)
Then the integral is evaluated and t is a scalar value. Mathematically, nothing has reallyhappened: g and f represent the same scalar value. Computationally, however, we have gonefrom a symbolic object which evaluates to a scalar, to the scalar itself.So, what happens in the following case? h = assemble(f_1) Users more frequently write a = ufl.TestFunction(V) , but that is simply syntactic sugar. A test function isjust argument 0 to a form. Equivalently, b = ufl.TrialFunction(W) . _1 contains the argument, or unknown Function , a . In other words, we need to evaluate h ( a ) = Z a d x = Z a i φ i d x (16)for unknown values a i . However the a i are scalar so they can be lifted out of the integral: Z a d x = Z φ i d x a i (17)Observe that we can now evaluate all the integrals, as they no longer contain unknown values.Next, we use the fact that φ ∗ i ( φ j ) = δ ij = I, the identity matrix. Hence: h ( a ) = Z φ i d x I ij a j = Z φ i d x φ ∗ i ( φ j ) a j = Z φ i d x φ ∗ i ( a j φ j )= Z φ i d x φ ∗ i ( a ) (18)We used the linearity of the dual basis functionals for the penultimate line. We have nowexpressed h in terms of the basis for V ∗ : h = h i φ ∗ i where: h i = Z φ i d x (19)The consequence of this is that the value of h is encoded in a vector of numbers, of the samesize as that which encodes c . Now one could say that the only purpose of assembly is to createa linear system to be solved, so the vector of values is the only thing which is required. Thisis the current FEniCS approach. This approach is problematic because it neglects the fact thatassembled operators can and are used by users in other contexts.What is missing, of course, is a reference to the dual basis { φ ∗ i } . The challenge is thatUFL does not have a class encoding { φ ∗ i } , and hence also does not have a suitable class for h .The Firedrake approach is to appeal to the natural isomorphism between V and V ∗ and simplymake h a Function . This has the advantage of preserving the link between the numerical valuesand the function space, but at the expense of incorrectly labelling t as being in V . What weshould really do is associate the list of coefficients h i with an object encoding the dual basis, ina manner analogous to the association between a list of coefficients and the basis for a functionspace in a Function .To illustrate the manner in which the assembled form encodes the same information as thesymbolic form, observe that: h ( c ) = f ( c ) = Z c d x = f = g (20)However the evaluation of h ( c ) is actually accomplished like this: h ( c ) = h i φ ∗ i ( c j φ j )= h i φ ∗ i ( φ j ) c j = h i δ ij c j = h i c i (21)So the assembly has turned the evaluation of an integral into the dot product between twovectors. 6 .2 UFL dual spaces and cofunctions The previous section suggests an outline of at least some missing UFL functionality. Thereshould be an object representing the dual to a function space (
DualFunctionSpace or maybejust
DualSpace ?) Presumably one would add a property to ufl.FunctionSpace to access thedual, and of course
DualSpace would also have this property but it would point back to theprimal space. Maybe ufl.FunctionSpace and ufl.DualSpace have a common parent class ufl.VectorSpace .We also need to be able to define a concrete member of V ∗ . Consistency with the existingUFL usage would suggest that this class would be ufl.Cocoefficient but this is a ridiculousand cumbersome name so perhaps ufl.Cofunction would be better. Firedrake and FEniCSsubclasses of ufl.Cofunction would, in a manner analogous to Function have both a vectorof values and a reference to the ufl.DualSpace . The question of which operations should bevalid for a cofunction requires more thought, though these would certainly include the basicvector operations of addition and multiplication by a scalar.
A cofunction is mathematically just a 1-form. The distinction between this and an existingUFL 1-form is that the latter is a symbolic expression that represents a 1-form, while the for-mer is a numerical 1-form represented directly in the dual basis. Since they are mathemati-cally objects of the same type, they should be arithmetically interchangeable. In particular, if f = ufl.Cofunction(V.dual) and v = ufl.Argument(V, 0) then the following makes math-ematical sense as a form and should be allowed: g = f + v*dx
This possibly implies that cofunctions intrinsically have an argument drawn from the relevantprimal space. Assembly of g can be achieved by simply defining assembly of a cofunction asthe identity operation. What do we get from:
A = assemble(f_2)
Mathematically this is defined as follows: A ij = Z φ i ψ j d x (22)That is to say, A is a matrix. Recall that { φ i } is the basis for V , which is the function space ofargument 0 in this case, and { ψ j } is the basis for W , which is argument 1. This means that theassembled operator has its arguments in ascending order, while the form notation has them indescending order. This is a frequent cause of confusion about which we will have to remainaware. FEniCS and Firedrake both agree that assembling a 2-form produces a matrix, thoughUFL does not have a corresponding class, so currently matrices are outside its scope. However,a matrix is an assembled 2-form in just the same way that a cofuncion is an assembled 1-form,so the logical conclusion would be the creation of a ufl.Matrix class which would behave asan assembled 2-form, and would have algebraic properties analogous to a cofunction. A matrixwould presumably have two intrinsic arguments.7 The delta operator and dual space arguments
Let’s now return to the interpolation operator. If we wish to incorporate it into UFL as a firstclass citizen then we need to take our extensions of the language just a little further. Let’ssuppose we can define arguments in dual spaces. These are mathematical objects of a dif-ferent kind to the arguments we have already met: they are unknown members of V ∗ ratherthan V . We suppose, therefore, that they are of a new type: ufl.Coargument . In order tominimise syntactic clutter, one might imagine overloading the object constructor such that ufl.Argument(V.dual, 0) is valid and returns a ufl.Coargument .This enables us to define a new form ufl.Delta such that if e is a UFL expression and v_ = ufl.Argument(V.dual, 0) then delta(e, v_) is a symbolic expression for the interpo-lation of e into V . If we assume that e does not contain any arguments and execute: f = assemble(Delta(e, v_)) Then f will be a Function in V representing the result of this interpolation. Conversely, if u = Argument (W, 1) then: A = assemble(Delta(u, v_)) results in a matrix which is the interpolation operator from W into V . One can also constructthe adjoint of this operator by reversing the argument positions: v = ufl.Argument(W, 0)u_ = ufl.Argument(V.dual, 1)AT = assemble(ufl.Delta(v, u_)) As a particular case of interest, if we write: u = ufl.Argument(V, 1)v_ = ufl.Argument(V.dual, 0)I = assemble(ufl.Delta(u, v_))
Then I is the identity matrix for V For completeness, we can also address the other cases. If we write: f = firedrake.Cofunction(V.dual)c = firedrake.Function(W)t = ufl.Delta(c, f)s = assemble(t) then t is a symbolic expression for f ( c ) and s is the scalar which results from evaluating it.Indeed, one might want to implement the __call__ special method on ufl.Cofunction and ufl.Coargument so that the syntax f(c) works directly. Note that this would simply return ufl.Delta(c, f) . ufl.Delta is still needed as the object which encodes the symbolic opera-tion of dual evaluation.Finally, if we write: f = firedrake.Cofunction(V.dual)u = ufl.Argument(W, 0)g = assemble(ufl.Delta(u, f)) Then g ∈ W ∗ is the interpolation of f into W ∗ . (Note that Cofunction interpolation is thetranspose of function interpolation). 8 Adjoint
The dual or adjoint of a 2-form is the conjugate of that form with its arguments relabelled inreversed order. This carries over to forms which map from or to dual spaces. If we take theexample of the delta form again, then: δ : V × W ∗ → R . (23)Equivalently: δ : V → W . (24)If we write δ ∗ for the adjoint to δ then: δ ∗ : W ∗ × V → R , (25)which is to say: δ ∗ : W ∗ → V ∗ . (26)Which is simply another way of saying that cofunction interpolation is the dual of coefficientinterpolation. The delta form is simply one example of an form whose first argument is in a dual space. Infact any operator which outputs a
Function and which is linear in any other arguments it has isa form of similar type. Note that such an operator could have known inputs (e.g. coefficients)with respect to which it is nonlinear, it is only the unknown arguments with respect to whichit has to be linear. In this context it is important to distinguish between the inputs to such anoperator, which we will call its operands and arguments. Arguments may appear in operands,but an operand also contain other UFL expressions. For example an operator N might takesome UFL expression as an input, evaluate a neural net and return the output as a Function .That is N ( e , m ) ∈ V for a UFL expression e , m ∈ R k and V some function space. Of coursesaying that N ( e , m ) → V is the same as saying N ( e , m ) × V ∗ → R , so an external operator isalso a form.
10 Form composition
Recall that a k -form of type V k − × . . . × V → R can be interpreted as an operator V k − × . . . × V → V ∗ . In other words, it can mathematically be used anywhere were an object in V ∗ isexpected. Similarly, the reflexivity of the vector spaces we work with means that a k -form oftype V k − × . . . × V ∗ → R can be interpreted as an operator V k − × . . . × V → V . That is, itcan be used where a value in V is expected.Hitherto, UFL forms have only taken inputs from primal spaces and produced outputs indual spaces, so composing forms by using the output of one as an input to another was not anavailable option. However, now that there are forms which produce primal space outputs, andforms which accept dual space inputs, we need to examine the consequences of this.Consider this simple case: v_ = ufl.Argument(V.dual, 0)f = firedrake.Function(W)v = ufl.Argument(V, 0)F_1 = ufl.Delta(f, v_) * v * dxl = assemble(F_1) W is defined on adifferent mesh from V . The assembly algorithm for l would be: tmp = firedrake.Function(V)l = assemble(ufl.replace(F_1, {ufl.Delta(f, w_): tmp})) I.e., we first evaluate the delta, and then substitute the result into the surrounding form. Thetype of the temporary variable can be inferred from the arguments to the delta. Consider nowthe 2-form case: v_ = ufl.Argument(V.dual, 0)w = ufl.Argument(W, 1)v = ufl.Argument(V, 0)F_2 = ufl.Delta(w, v_) * v * dxA = assemble(F_2)
This results in the following assembly algorithm: tmp_1 = assemble(Delta(w, v_)) tmp_2 = assemble(ufl.replace(F_2, {ufl.Delta(w, v_): ufl.Argument(V, 1)}))A = tmp_2 @ tmp_1
Again the, the types of the temporary variables can be inferred from the arguments of the delta.The delta has 2 arguments and they are in the canonical order, so it is replaced by an argumentnumbered 1 but in the function space dual to the argument 0 function space. The sequence ofthe matrix multiplication is similarly inferable from the arguments and the sequence of com-position (the delta occurs inside the other form, and not vice versa). The general rule is thatthe form is replaced by an argument in the dual space to its zero argument, and with an argu-ment number which is the next available argument number in the form into which it is beingsubstituted.
What are the arguments of
F_2 ? One might be tempted to say v_ , w , and v , but this wouldcharacterise F_2 as a 3-form rather than a 2-form. It would also have two different argumentsnumbered 0 ( v_ , and v ), so it would not even be a well-formed form. The answer to this puzzleis that the 0 argument of a form nested into another form is consumed by the substitution. Inthis case, this means that the arguments of the composite form are v and w . Consider instead the case where the interpolation happens for the test function. By symmetrywith
F_2 , one might expect to be able to write: w_ = ufl.Argument(W.dual, 0)v = ufl.Argument(V, 0)u = ufl.Argument(V, 1)F_3 = u * ufl.Delta(v, w_) * dx
However this is ill-formed since the delta contains two arguments numbered 0. This difficultyis the result of the fact that form composition only makes sense if we think of forms as operatorsinto the dual space of argument 0. In this context, argument 0 is not a form input at all, but isinstead a placeholder for the result of evaluating the form for known values of all of its higher-numbered arguments. However, composition of linear operators works in both directions sowe can instead write: 10 = ufl.Argument(W, 0)v = ufl.Argument(W, 0)u = ufl.Argument(V, 1)F_4 = ufl.Delta(w, u * v * dx)
Which is well-formed UFL, and mathematically equivalent to what
F_3 attempted to achieve.As one might expect, F : V → W ∗ .
11 Action
The action of a form on a coefficient or argument simply replaces the highest numbered argu-ment with that coefficient. This carries over directly to the dual case. If the highest numberedargument of a form lies in a dual space, then the action is only well defined on a cofunction orcoargument from the appropriate space.Form composition means that it should also be possible to take the action of a form onanother form, so long as the space of the highest numbered argument of the first form is dualto the space of the 0 argument to the second form. That is, an alternative definition of
F_4 isgiven by: w = ufl.Argument(W, 0)omega = ufl.Argument(W.dual, 1)v = ufl.Argument(V, 0)u = ufl.Argument(V, 1)F_4 = ufl.action(ufl.Delta(w, omega), u * v * dx)
Of course this lines up with the assembly strategy for
F_4 , which would be: tmp_1 = assemble(Delta(w_, ufl.Argument(V.dual, 1))tmp_2 = assemble(u * v * dx)A = tmp_1 @ tmp_2
References
Alnæs, Martin S et al. (2014). “Unified form language: A domain-specific language for weakformulations of partial differential equations”. In:
ACM Transactions on Mathematical Soft-ware (TOMS)