%matplotlib inline
import matplotlib.pyplot as plt
import firedrake
import icepack, icepack.plot

Meshes, functions

This demo will show some of the basic features of the finite element modeling package firedrake, which icepack is built on. These are the kinds of things you’d need to know before you begin to model anything. First, we’ll show to create a simple mesh, and how to define a function on that mesh symbolically using handles for the x- and y-coordinates of the points of the domain. Then we’ll demonstrate some of the plotting routines included in icepack, which are meant to mimic as closely as possible the corresponding routines in matplotlib. Finally, we’ll show how to analyze a function by evaluating some area and contour integrals.


First, we need to make a mesh. Firedrake includes routines for generating meshes of simple domains; to see all of the meshes that are available, you can type


from the python interpreter. In this demo we’ll use a mesh of the unit square.

nx, ny = 16, 16
mesh = firedrake.UnitSquareMesh(nx, ny)

Mesh plots are shown with a different color for each segments of the boundary. You can then create a legend that will show you how the colors correspond to the numeric IDs of each boundary segment. This is useful to when applying different boundary conditions to different segments, for example, Dirichlet boundary conditions on part of the domain and Neumann boundary conditions on the other part.

fig, axes = icepack.plot.subplots()
icepack.plot.triplot(mesh, axes=axes)


We can make scalar and vector fields defined on this mesh by interpolating the values of an algebraic expression to some space of functions.

First, we have to construct a function space Q. We already have a mesh, and we need to decide what element family and polynomial degree. In almost all the cases you’ll encounter, the element family consists of continuous piecewise-polynomial functions in each triangle, which is abbreviated to CG for “continuous Galerkin”.

Next, we have to make an expression for the function and interpolate it to the function space Q. The function firedrake.SpatialCoordinate returns two symbolic objects x, y for the coordinates of each point of the mesh. We can then use these symbols to define expressions for the function we’d like to analyze. In this case, I’ve chosen the Rosenbrock function. Firedrake has built-in functions for evaluating various transcendental functions of the coordinates, for example the sine, cosine, exponential, logarithm, etc. To see all of the available functions, you can check the namespace ufl.operators.

Finally, the function firedrake.interpolate takes in an expression and a function space, and returns a field from that function space.

Q = firedrake.FunctionSpace(mesh, family='CG', degree=2)

x, y = firedrake.SpatialCoordinate(mesh)
a, b = 0.5, 10.0
expr = (a - x)**2 + b*(y - x**2)**2
q = firedrake.interpolate(expr, Q)

fig, axes = icepack.plot.subplots()
contours = icepack.plot.tricontourf(q, 36, axes=axes)

We can make vector fields in much the same way as we make scalar fields, the only difference being the type of function space. There are many ways to plot a vector field. Here we show a plot of the streamlines of the vector field, colored according to the magnitude of the vector. This method is expensive, especially if you set the resolution or spacing of the streamlines to be too small, but nonetheless it produces really nice graphics.

V = firedrake.VectorFunctionSpace(mesh, family='CG', degree=2)
v = firedrake.interpolate(firedrake.grad(expr), V)

fig, axes = icepack.plot.subplots()
opts = {'precision': 1/128, 'density': 1/16, 'max_num_points': 128}
streamlines = icepack.plot.streamplot(v, axes=axes, **opts)


We’ve already seen about how to visualize a scalar or a vector field; this is just a call to the function icepack.plot. To analyze the results of a simulation, you might also want to evaluate scalar or vector fields at points in the domain:

z = (a, a**2)
[ 3.46901282e-15 -2.52175273e-15]

Firedrake also provides a rich set of operations for evaluating integral expressions of scalar and vector fields. Much like how the function SpatialCoordinate gives you two symbols x, y that represent the coordinates of each point of the mesh, firedrake also provides an object dx that represents the differential area element.

To define an integral, we multiply an expression by dx, and then call the function firedrake.assemble to evaluate it.

from firedrake import inner, grad, dx, assemble

print(assemble(inner(v, v)*dx))

We can also evaluate integrals over the boundary of the mesh. Boundary integrals are specified using the differential surface element ds instead of dx.

Sometimes we need to integrate over only part of the boundary. We can specify which part by passing the right numeric ID to ds; these are the same IDs that are color-coded in the mesh plot above.

Finally, you might also want to evaluate the flux of a vector field across the boundary. To evaluate a flux, we need to get the unit outward normal vector to the boundary. You can get a symbol representing the unit outward normal by calling the function firedrake.FacetNormal.

from firedrake import ds


n = firedrake.FacetNormal(mesh)
print(assemble(inner(v, n) * ds))


Icepack uses the functions in firedrake to implement solvers for the various physics problems that show up in ice sheet modeling. The firedrake routines that we’ve demonstrated here can be used for analyzing simulation results, either through visualization, evaluating fields at points, or evaluating integrals of fields. They can also be used to define input fields of simulations, if those fields have a simple analytical expression. In the next demo, we’ll show how to use icepack to solve for the velocity and thickness of a floating ice shelf with a synthetic geometry. In later demos, we’ll show how to use real observational data sets, as well as more complicated geometries imported from mesh generators.

To read more about firedrake, you can visit their documentation or check out some of the demos.