Contents

# Visualising the results of simulations¶

Having run a simulation, it is likely that we will want to look at the results. To do this, Firedrake supports saving data in VTK format, suitable for visualisation in Paraview (amongst others).

In addition, 1D and 2D function could be plotted and displayed using the python library of matplotlib (an optional dependency of firedrake)

## Creating output files¶

Output for visualisation purposes is managed with a `File`

object. To create one, we just need to pass the name of the output
file on disk. The file Firedrake creates is in PVD and therefore the
requested file name must end in `.pvd`

.

```
outfile = File("output.pvd")
# The following raises an error
badfile = File("output.vtu")
```

To save functions to the `File`

we use the
`write()`

method.

```
mesh = UnitSquareMesh(1, 1)
V = FunctionSpace(mesh, "DG", 0)
f = Function(V)
f.interpolate(sin(SpatialCoordinate(mesh)[0]))
outfile = File("output.pvd")
outfile.write(f)
```

Note

Output created for visualisation purposes is not intended for purposes other than visualisation. If you need to save data for checkpointing purposes, you should instead use Firedrake’s checkingpointing capabilities.

## Saving time-dependent data¶

Often, we have a time-dependent simulation and would like to save the
same function at multiple timesteps. This is straightforward, we must
create the output `File`

outside the time loop and call
`write()`

inside.

```
...
outfile = File("timesteps.pvd")
while t < T:
...
outfile.write(f)
t += dt
```

The PVD data format supports specifying the timestep value for
time-dependent data. We do not have to provide it to
`write()`

, by default an integer counter is used that is
incremented by 1 each time `write()`

is called. It is
possible to override this by passing the keyword argument `time`

.

```
...
outfile = File("timesteps.pvd")
while t < T:
...
outfile.write(f, time=t)
t += dt
```

## Visualising high-order data¶

The file format Firedrake outputs to currently supports the
visualisation of scalar-, vector-, or tensor-valued fields
represented with an arbitrary order (possibly discontinuous) Lagrange basis.
Furthermore, the fields must be in an isoparametric function space, meaning
the mesh coordinates associated to a field must be represented
with the same basis as the field. To visualise fields in anything
other than these spaces we must transform the data to this
format first. One option is to do so by hand before outputting.
Either by interpolating or else `projecting`

the mesh coordinates and then the field. Since this is
such a common operation, the `File`

object is set up to manage these
operations automatically, we just need to choose whether we want data to be
interpolated or projected. The default is to use interpolation. For example,
assume we wish to output a vector-valued function that lives in an \(H(\operatorname{div})\)
space. If we want it to be interpolated in the output file we can use

```
V = FunctionSpace(mesh, "RT", 2)
f = Function(V)
...
outfile = File("output.pvd")
outfile.write(f)
```

If instead we want projection, we use

```
projected = File("proj_output.pvd", project_output=True)
projected.write(f)
```

Note

This feature requires Paraview version 5.5.0 or better. If you must use an
older version of Paraview, you must manually interpolate mesh coordinates
and field coordinates to a piecewise linear function space, represented
with either a Lagrange (H1) or discontinuous Lagrange (L2) basis. The `File`

is also setup to manage this issue. For instance, we can force the output
to be discontinuous piecewise linears via

```
projected = File("proj_output.pvd", target_degree=1, target_continuity=H1)
projected.write(f)
```

### Using Paraview on higher order data¶

Paraview’s visualisation algorithims are typically exact on piecewise linear data, but if you write higher order data, Paraview will produce an approximate visualisation. This approximation can be controlled in at least two ways:

Under the display properties of an unstructured grid, the Nonlinear Subdivision Level can be increased; this option controls the display of unstructured grid data and can be used to present a plausible curved geometry. Further, the Nonlinear Subdivision Level can also be changed after applying filters such as Extract Surface.

The Tessellate filter can be applied to unstructured grid data and has two parameters: Chord Error, Maximum Number of Subdivisions, and Field Error. Tessellation is the process of approximating a higher order geometry via subdividing cells into smaller linear cells. Chord Error is a tessellation error metric, the distance between the midpoint of any edge on the tessellated geometry and a corresponding point in the original geometry. Field Error is analogous to Chord Error: the error of the field on the tessellated data is compared pointwise to the original data at the midpoints of the edges of the tessellated geometry and the corresponding points on the original geometry. The Maximum Number of Subdivisions is the maximum number of times an edge in the original geometry can be subdivided.

Besides the two tools listed above, Paraview provides many other tools (filters) that might be applied to the original data or composed with the tools listed above. Documentation on these interactions is sparse, but tessellation can be used to understand this issue: the Tessellate filter produces another unstructured grid from its inputs so algorithms can be applied to both the tessellated and input unstructured grid. The tessellated data can also be saved for future reference.

Note

Field Error is hidden in the current Paraview UI (5.7) so we include a visual guide wherein the field error is set via the highlighted field directly below Chord Error:

We also note that the Tessellate filter (and other filters) can be more clearly controlled via the Paraview Python shell (under the View menu). For instance, Field Error can be more clearly specified via an argument to the Tessellate filter constructor.

```
from paraview.simple import *
pvd = PVDReader(FileName="Example.pvd")
tes = Tessellate(pvd, FieldError=0.001)
```

## Saving multiple functions¶

Often we will want to save, and subsequently visualise, multiple
different fields from a simulation. For example the velocity and
pressure in a fluids models. This is possible either by having a
separate output file for each field, or by saving multiple fields to
the same output file. The latter may be more convenient for
subsequent analysis. To do this, we just need to pass multiple
`Function`

s to `write()`

.

```
u = Function(V, name="Velocity")
p = Function(P, name="Pressure")
outfile = File("output.pvd")
outfile.write(u, p, time=0)
# We can happily do this in a timeloop as well.
while t < t:
...
outfile.write(u, p, time=t)
```

Note

Subsequent writes to the same file *must* use the same number of
functions, and the functions must have the *same* names. The
following example results in an error.

```
u = Function(V, name="Velocity")
p = Function(P, name="Pressure")
outfile = File("output.pvd")
outfile.write(u, p, time=0)
...
# This raises an error
outfile.write(u, time=1)
# as does this
outfile.write(p, u, time=1)
```

### Selecting the output space when outputting multiple functions¶

All functions, including the mesh coordinates, that are output to the same file must be represented in the same space, the rules for selecting the output space are as follows. First, all functions must be defined via the same cell type otherwise an exception will be thrown. Second, if all functions are continuous (i.e. they live in \(H^1\)), then the output space will be a piecewise continuous space. If any of the functions are at least partially discontinuous, again including the coordinate field (this occurs when using periodic meshes), then the output will use a piecewise discontinuous space. Third, the degree of the basis will be the maximum degree used over the spaces of all input functions. For elements where the degree is a tuple (this occurs when using tensor product elements), the the maximum will be over the elements of the tuple too, meaning a tensor product of elements of degree 4 and 2 will be turned into a tensor product of elements of degree 4 and 4.

## Plotting with matplotlib¶

Plotting 1D and 2D functions could be as easy as calling the built-in plot
function `plot`

with the `Function`

you wish to
plot.

Currently, firedrake supports plotting 1D and 2D functions, this is made possible with an optional dependency matplotlib package.

To install matplotlib, please look at the installation instructions of matplotlib.

For 1D functions with degree less than 4, the plot of the function would be
exact using Bezier curves. For higher order 1D functions, the plot would be the
linear approximation by sampling points of the function. The number of sample
points per element could be specfied to when calling `plot`

.

For multiple 1D functions, for example, in the case of time-dependent functions
at different times. They could be plotted together by passing the list of
function when calling the function `plot`

. The returned
figure will contain a slider and an autoplay button so that it could be viewed
in a animated fashion. The plus and minus buttons can change the speed of the
animation.

When used in Jupyter Notebook, plotting multiple 1D functions using additional
keyword argument `interactive=True`

when calling the function
`plot`

will generate an interactive slider for
selecting the figures.

For 2D functions, both surface plots and contour plots are supported. By
default, the `plot`

will return a surface plot in the
colour map of coolwarm. Contour plotting could be enabled by passing the keyword
argument `contour=True`

.