Algorithms

Algorithms to work with UFL forms and expressions can be found in the submodule ufl.algorithms. You can import all of them with the line

from ufl.algorithms import *

This chapter gives an overview of (most of) the implemented algorithms. The intended audience is primarily developers, but advanced users may find information here useful for debugging.

While domain specific languages introduce notation to express particular ideas more easily, which can reduce the probability of bugs in user code, they also add yet another layer of abstraction which can make debugging more difficult when the need arises. Many of the utilities described here can be useful in that regard.

Formatting expressions

Expressions can be formatted in various ways for inspection, which is particularly useful for debugging. We use the following as an example form for the formatting sections below:

element = FiniteElement("CG", triangle, 1)
v = TestFunction(element)
u = TrialFunction(element)
c = Coefficient(element)
f = Coefficient(element)
a = c*u*v*dx + f*v*ds

str

Compact, human readable pretty printing. Useful in interactive Python sessions. Example output of str(a):

{ v_0 * v_1 * w_0 } * dx(<Mesh #-1 with coordinates parameterized by <Lagrange vector element of degree 1 on a triangle: 2 x <CG1 on a triangle>>>[everywhere], {})
+  { v_0 * w_1 } * ds(<Mesh #-1 with coordinates parameterized by <Lagrange vector element of degree 1 on a triangle: 2 x <CG1 on a triangle>>>[everywhere], {})

repr

Accurate description of an expression, with the property that eval(repr(a)) == a. Useful to see which representation types occur in an expression, especially if str(a) is ambiguous. Example output of repr(a):

Form([Integral(Product(Argument(FunctionSpace(Mesh(VectorElement('Lagrange', triangle, 1, dim=2), -1), FiniteElement('Lagrange', triangle, 1)), 0, None), Product(Argument(FunctionSpace(Mesh(VectorElement('Lagrange', triangle, 1, dim=2), -1), FiniteElement('Lagrange', triangle, 1)), 1, None), Coefficient(FunctionSpace(Mesh(VectorElement('Lagrange', triangle, 1, dim=2), -1), FiniteElement('Lagrange', triangle, 1)), 0))), 'cell', Mesh(VectorElement('Lagrange', triangle, 1, dim=2), -1), 'everywhere', {}, None), Integral(Product(Argument(FunctionSpace(Mesh(VectorElement('Lagrange', triangle, 1, dim=2), -1), FiniteElement('Lagrange', triangle, 1)), 0, None), Coefficient(FunctionSpace(Mesh(VectorElement('Lagrange', triangle, 1, dim=2), -1), FiniteElement('Lagrange', triangle, 1)), 1)), 'exterior_facet', Mesh(VectorElement('Lagrange', triangle, 1, dim=2), -1), 'everywhere', {}, None)])

Tree formatting

ASCII tree formatting, useful to inspect the tree structure of an expression in interactive Python sessions. Example output of tree_format(a):

Form:
  Integral:
    integral type: cell
    subdomain id: everywhere
    integrand:
      Product
      (
        Argument(FunctionSpace(Mesh(VectorElement('Lagrange', triangle, 1, dim=2), -1), FiniteElement('Lagrange', triangle, 1)), 0, None)
        Product
        (
          Argument(FunctionSpace(Mesh(VectorElement('Lagrange', triangle, 1, dim=2), -1), FiniteElement('Lagrange', triangle, 1)), 1, None)
          Coefficient(FunctionSpace(Mesh(VectorElement('Lagrange', triangle, 1, dim=2), -1), FiniteElement('Lagrange', triangle, 1)), 0)
        )
      )
Integral:
  integral type: exterior_facet
  subdomain id: everywhere
  integrand:
    Product
    (
      Argument(FunctionSpace(Mesh(VectorElement('Lagrange', triangle, 1, dim=2), -1), FiniteElement('Lagrange', triangle, 1)), 0, None)
      Coefficient(FunctionSpace(Mesh(VectorElement('Lagrange', triangle, 1, dim=2), -1), FiniteElement('Lagrange', triangle, 1)), 1)
    )

Inspecting and manipulating the expression tree

This subsection is mostly for form compiler developers and technically interested users.

Traversing expressions

iter_expressions

Example usage:

for e in iter_expressions(a):
    print str(e)

outputs:

v_0 * v_1 * w_0
v_0 * w_1

Transforming expressions

So far we presented algorithms meant to inspect expressions in various ways. Some recurring patterns occur when writing algorithms to modify expressions, either to apply mathematical transformations or to change their representation. Usually, different expression node types need different treatment.

To assist in such algorithms, UFL provides the Transformer class. This implements a variant of the Visitor pattern to enable easy definition of transformation rules for the types you wish to handle.

Shown here is maybe the simplest transformer possible:

class Printer(Transformer):
    def __init__(self):
        Transformer.__init__(self)

    def expr(self, o, *operands):
        print "Visiting", str(o), "with operands:"
        print ", ".join(map(str,operands))
        return o

element = FiniteElement("CG", triangle, 1)
v = TestFunction(element)
u = TrialFunction(element)
a = u*v

p = Printer()
p.visit(a)

The call to visit will traverse a and call Printer.expr on all expression nodes in post–order, with the argument operands holding the return values from visits to the operands of o. The output is:

Visiting v_0 * v_1 with operands:
v_0, v_1

\((v^0_h)(v^1_h)\)

Implementing expr above provides a default handler for any expression node type. For each subclass of Expr you can define a handler function to override the default by using the name of the type in underscore notation, e.g. vector_constant for VectorConstant. The constructor of Transformer and implementation of Transformer.visit handles the mapping from type to handler function automatically.

Here is a simple example to show how to override default behaviour:

from ufl.classes import *
class CoefficientReplacer(Transformer):
    def __init__(self):
        Transformer.__init__(self)

    expr = Transformer.reuse_if_possible
    terminal = Transformer.always_reuse

    def coefficient(self, o):
        return FloatValue(3.14)

element = FiniteElement("CG", triangle, 1)
v = TestFunction(element)
f = Coefficient(element)
a = f*v

r = CoefficientReplacer()
b = r.visit(a)
print b

which outputs

3.14 * v_0

The output of this code is the transformed expression b == 3.14*v. This code also demonstrates how to reuse existing handlers. The handler Transformer.reuse_if_possible will return the input object if the operands have not changed, and otherwise reconstruct a new instance of the same type but with the new transformed operands. The handler Transformer.always_reuse always reuses the instance without recursing into its children, usually applied to terminals. To set these defaults with less code, inherit ReuseTransformer instead of Transformer. This ensures that the parts of the expression tree that are not changed by the transformation algorithms will always reuse the same instances.

We have already mentioned the difference between pre–traversal and post–traversal, and some times you need to combine the two. Transformer makes this easy by checking the number of arguments to your handler functions to see if they take transformed operands as input or not. If a handler function does not take more than a single argument in addition to self, its children are not visited automatically, and the handler function must call visit on its operands itself.

Here is an example of mixing pre- and post-traversal:

class Traverser(ReuseTransformer):
    def __init__(self):
        ReuseTransformer.__init__(self)

    def sum(self, o):
        operands = o.operands()
        newoperands = []
        for e in operands:
            newoperands.append( self.visit(e) )
        return sum(newoperands)

element = FiniteElement("CG", triangle, 1)
f = Coefficient(element)
g = Coefficient(element)
h = Coefficient(element)
a = f+g+h

r = Traverser()
b = r.visit(a)
print b

This code inherits the ReuseTransformer as explained above, so the default behaviour is to recurse into children first and then call Transformer.reuse_if_possible to reuse or reconstruct each expression node. Since sum only takes self and the expression node instance o as arguments, its children are not visited automatically, and sum explicitly calls self.visit to do this.

Automatic differentiation implementation

This subsection is mostly for form compiler developers and technically interested users.

First of all, we give a brief explanation of the algorithm. Recall that a Coefficient represents a sum of unknown coefficients multiplied with unknown basis functions in some finite element space.

\[w(x) = \sum_k w_k \phi_k(x)\]

Also recall that an Argument represents any (unknown) basis function in some finite element space.

\[v(x) = \phi_k(x), \qquad \phi_k \in V_h .\]

A form \(L(v; w)\) implemented in UFL is intended for discretization like

\[b_i = L(\phi_i; \sum_k w_k \phi_k), \qquad \forall \phi_i \in V_h .\]

The Jacobi matrix \(A_{ij}\) of this vector can be obtained by differentiation of \(b_i\) w.r.t. \(w_j\), which can be written

\[A_{ij} = \frac{d b_i}{d w_j} = a(\phi_i, \phi_j; \sum_k w_k \phi_k), \qquad \forall \phi_i \in V_h, \quad \forall \phi_j \in V_h ,\]

for some form a. In UFL, the form a can be obtained by differentiating L. To manage this, we note that as long as the domain \(\Omega\) is independent of \(w_j\), \(\int_\Omega\) commutes with \(\frac{d}{d w_j}\), and we can differentiate the integrand expression instead, e.g.,

\[\begin{split}L(v; w) = \int_\Omega I_c(v; w) \, dx + \int_{\partial\Omega} I_e(v; w) \, ds, \\ \frac{d}{d w_j} L(v; w) = \int_\Omega \frac{d I_c}{d w_j} \, dx + \int_{\partial\Omega} \frac{d I_e}{d w_j} \, ds.\end{split}\]

In addition, we need that

\[\frac{d w}{d w_j} = \phi_j, \qquad \forall \phi_j \in V_h ,\]

which in UFL can be represented as

\[\begin{split}w &= \mathtt{Coefficient(element)}, \\ v &= \mathtt{Argument(element)}, \\ \frac{dw}{d w_j} &= v,\end{split}\]

since \(w\) represents the sum and \(v\) represents any and all basis functions in \(V_h\).

Other operators have well defined derivatives, and by repeatedly applying the chain rule we can differentiate the integrand automatically.

Computational graphs

This section is for form compiler developers and is probably of no interest to end-users.

An expression tree can be seen as a directed acyclic graph (DAG). To aid in the implementation of form compilers, UFL includes tools to build a linearized [1] computational graph from the abstract expression tree.

A graph can be partitioned into subgraphs based on dependencies of subexpressions, such that a quadrature based compiler can easily place subexpressions inside the right sets of loops.

[1]

Linearized as in a linear datastructure, do not confuse this with automatic differentiation.

The computational graph

Consider the expression

\[f = (a + b) * (c + d)\]

where a, b, c, d are arbitrary scalar expressions. The expression tree for f looks like this:

a   b   c   d
\   /   \   /
  +       +
   \     /
      *

In UFL f is represented like this expression tree. If a, b, c, d are all distinct Coefficient instances, the UFL representation will look like this:

Coefficient Coefficient Coefficient Coefficient
          \     /             \     /
            Sum                 Sum
               \               /
                --- Product ---

If we instead have the expression

\[f = (a + b) * (a - b)\]

the tree will in fact look like this, with the functions a and b only represented once:

Coefficient     Coefficient
   |       \   /       |
   |        Sum      Product -- IntValue(-1)
   |         |         |
   |       Product     |
   |         |         |
   |------- Sum -------|

The expression tree is a directed acyclic graph (DAG) where the vertices are Expr instances and each edge represents a direct dependency between two vertices, i.e. that one vertex is among the operands of another. A graph can also be represented in a linearized data structure, consisting of an array of vertices and an array of edges. This representation is convenient for many algorithms. An example to illustrate this graph representation follows:

G = V, E
V = [a, b, a+b, c, d, c+d, (a+b)*(c+d)]
E = [(6,2), (6,5), (5,3), (5,4), (2,0), (2,1)]

In the following, this representation of an expression will be called the computational graph. To construct this graph from a UFL expression, simply do

G = Graph(expression)
V, E = G

The Graph class can build some useful data structures for use in algorithms:

Vin  = G.Vin()  # Vin[i]  = list of vertex indices j such that there is an edge from V[j] to V[i]
Vout = G.Vout() # Vout[i] = list of vertex indices j such that there is an edge from V[i] to V[j]
Ein  = G.Ein()  # Ein[i]  = list of edge indices j such that E[j] is an edge to V[i], e.g. E[j][1] == i
Eout = G.Eout() # Eout[i] = list of edge indices j such that E[j] is an edge from V[i], e.g. E[j][0] == i

The ordering of the vertices in the graph can in principle be arbitrary, but here they are ordered such that

\[v_i \prec v_j, \quad \forall j > i,\]

where \(a \prec b\) means that \(a\) does not depend on \(b\) directly or indirectly.

Another property of the computational graph built by UFL is that no identical expression is assigned to more than one vertex. This is achieved efficiently by inserting expressions in a dict (a hash map) during graph building.

In principle, correct code can be generated for an expression from its computational graph simply by iterating over the vertices and generating code for each one separately. However, we can do better than that.

Partitioning the graph

To help generate better code efficiently, we can partition vertices by their dependencies, which allows us to, e.g., place expressions outside the quadrature loop if they don’t depend (directly or indirectly) on the spatial coordinates. This is done simply by

P = partition(G)