Unstructured Multigrid

In a previous document, we derived an algorithm for coarsening simplex meshes, motivated by how it would enable us to apply multigrid methods to complicated geometries:


In this document, we'll look at how to use a hierarchy of coarsened meshes to implement some standard multigrid methods and see how well they work.

As a matter of notation, we'll refer to quantities on original mesh as h, where is some quantity, and h represents the element size. The analogous quantities on the coarser and coarser meshes will be referred to as 2h,4h,8h, and so on. Technically, the element size difference between our unstructured coarse meshes isn't exactly a factor of two like it is for structured multigrid, but we'll maintain this notation to be consistent with existing multigrid conventions.

Interpolation and Restriction

To start, we need a way to move information between meshes, so let's define some symbols to represent the important operations.

h:=I2hInterpolation operator2h2h:=R2hRestriction operatorh

Interpolation (also sometimes referred to as "prolongation" in multigrid) takes coarse quantities and moves them to a finer discretization while restriction does the opposite. Since there are fewer degrees of freedom on the coarse meshes, the matrices for the interpolation operators I will be "tall and skinny", and R will be "short and fat".

The interpolation operators can be defined unambigously for continuous finite element fields by recognizing that the mesh nodes themselves should satisfy


where X is a matrix whose rows are the nodes of a mesh. Finite element interpolation is done element-by-element, so we first have to figure out which coarse element contains a given node of the fine mesh. Then, we determine the isoparametric coordinates (ξ,η) for that point. An illustration of this is given below, where the coarse mesh is shown in black and the red point represents one arbitrary node from the fine mesh.


The isoparametric coordinates are used to evaluate that element's shape functions, which become the values of the nonzero entries in that row of I2h.

But what if some of the fine mesh nodes are not inside any element of the coarse mesh?

This edge case doesn't show up with cartesian structured meshes, but it can happen with our coarsened unstructured meshes. We follow the same procedure, except that instead of using the element that contains the point, use the element closest to it.


Finally, there are a few different ways to define restriction operators, but the one we'll use is the "variationally correct" choice of R2h:=(I2h). In addition to being simple, it also has the benefit of preserving symmetry (if present) of the coarse operators.

Mathematica Implementation

Note: this FindNearestTriangle[] function is a naive O(n2) implementation. For large meshes, you would want to use something like a BVH to accelerate this sort of spatial query.

Physics on the Coarse Meshes

The last thing we need is a way to represent the physics on the coarse meshes. There are two natural ways to do this:

  1. Use the interpolation/restriction operators to form A2h:=R2hAhI2h

    Accuracy: this choice has ideal spectral properties

    Memory Use: requires forming A2h,A4h,... explicitly

    Complexity: sparse matrix assembly is expensive and complicated

  2. Directly using the coarse meshes to discretize the problem

    Simple: reuses existing finite element tools

    Fast: allows use of optimized residual evaluation kernels

    Memory use: doesn't need to form A2h,A4h,... explicitly

    Accuracy: lacks the ideal spectral properties of option 1

Option 1 seems to be favored by the multigrid community, but option 2 has some important advantages from a software development perspective that make it appealing. Ultimately, the best choice is determined by how well they perform, so we need to do some numerical experiments to weigh the benefits and drawbacks of each approach. We'll use a Poisson problem defined on the wrench meshes shown at the top of the document, with original coefficient matrix given by


for the examples to follow. Righthand sides are randomly generated with uniformly random entries on the interval (1,1).

Multigrid Cycles

Now that we have our interpolation, restriction and physics defined over the entire hierarchy, we can finally start to implement some multigrid cycles. The two we'll consider in this document are represented pictorially below


In these diagrams, the height of each node represents its level in the mesh hierarchy (lower = coarser). The diagrams are read left-to-right, and each node along the path implies some calculation taking place at that level (typically smoothing). Lines between nodes mean moving data from one discretization to another, so it will involve either interpolation (if line's slope is positive) or restriction (if the slope is negative).

We can see that the V cycle starts at the finest discretization, and alternates smoothing and restriction until it reaches the bottom. At this coarsest level, we can either factorize the coarsest system and solve, or just apply smoothing. Following that, it interpolates and smooths at each level until coming back to the finest mesh.

In contrast, the full multigrid (FMG) cycle starts at the bottom of the hierarchy, and alternates interpolation and V cycles of increasing depth until reaching the finest mesh. Although the FMG cycle looks considerably "wider" than the V cycle of the same depth, the actual time to complete it is not considerably longer, as most of the extra steps occur at the coarsest levels, where calculations are inexpensive (especially in 3D).

Mathematica Implementation

These two methods can be implemented succinctly as recursive functions:

As Fixed-Point iterations

Fixed point iteration was the original iterative way of solving Ax=b. Starting with some x0, subsequent iterates are determined by


for some choice of matrix M (for example, Jacobi method takes M=diag(A)1). This method is beautiful in its simplicity, but it just doesn't work very well by itself. If we try applying Jacobi iterations to our system of equations and plot the relative residuals


we see that the residual drops quickly for the first few iterates, but convergence slows down and stagnates. Compare that to "M" corresponding to a 5-deep V cycle with k Jacobi iterations per level, and a direct solve at the coarsest level (44 degrees of freedom)


This is clearly a huge improvement over regular Jacobi method. Performing more smoothing at each stage in the V cycle makes each iteration more effective, but also more expensive. A V cycle with k = 4, may take as long as 10-15 Jacobi iterations, but it does a much better job decreasing the residual.

Let's also compare V cycles to FMG cycles, and look at the convergence of options 1 and 2 for the coarse representations of A2h,A4h,...


As we might expect, FMG cycles converges slightly faster than V cycles per iteration, but they also do more calculation. With the V cycles, option 2 lags behind slightly , taking about 20-30% more iterations to reach the same relative residual as option 1. However, there is practically zero distinction between options 1 and 2 with FMG.

As Preconditioners

Krylov methods are the more modern way of solving systems of equations iteratively. They use the subspace spanned by {b,Ab,A2b,A3b...Anb} to form approximate solutions to the system of equations. Much like fixed point iterations, naively applying a Krylov method often gives poor results, if we apply conjugate gradient (cg) to our system of equations (κ30,000,000):


it takes hundreds of iterations before making any real progress. Eventually, cg overtakes Jacobi, before stagnating for another few hundred cycles. Fortunately, Krylov methods can perform much more effectively and reliably with an appropriate choice of preconditioner. Instead of using M in a fixed point iteration, we can use it as a preconditioner instead


A surprisingly large number of real applications actually use (as a coworker puts it) "Jacobi and pray" as a preconditioning strategy, but it frequently performs poorly. The V cycle as a fixed point iteration performs quite well, but the multigrid-as-a-preconditioner options work so effectively that they can hardly even be seen on this plot. If we zoom in significantly we can resolve the differences


In just a handful of iterations, Krylov + multigrid-as-a-preconditioner options were able to drive the residual to machine precision.

Elasticity Example

This document is getting kind of long, so we'll continue the discussion of multigrid applied to elasticity on a separate page.


For this example problem, some key takeaways were that:

But this wasn't a particularly challenging problem, so it would be inappropriate to try and claim that these conclusions were true in general. That being said, the success of unstructured multigrid on this example is encouraging, and (in my mind) warrants more numerical experiments. Some of the questions I'd like to investigate include