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$$\square^h$, where $◻$$\square$ is some quantity, and $h$$h$ represents the element size. The analogous quantities on the coarser and coarser meshes will be referred to as ${◻}^{2h},{◻}^{4h},{◻}^{8h},$$\square^{2h}, \square^{4h}, \square^{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.

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 $\mathbf{I}$$\bold{I}$ will be "tall and skinny", and $\mathbf{R}$$\bold{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 $\mathbf{X}$$\bold{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 $(\xi ,\eta )$$(\xi, \eta)$ 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 ${\mathbf{I}}^{2h}$$\bold{I}^{2h}$.

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 ${\mathbf{R}}^{2h}:=({\mathbf{I}}^{2h}{)}^{\mathrm{\top }}$$\bold{R}^{2h} := (\bold{I}^{2h})^\top$. In addition to being simple, it also has the benefit of preserving symmetry (if present) of the coarse operators.

Mathematica Implementation

IsoparametricCoordinates[p_, tri_] := Module[{A, b},
  A = Transpose[{tri[[2]] - tri[[1]], tri[[3]] - tri[[1]]}];
  b = p - tri[[1]];
  LinearSolve[A, b]
]

FindNearestTriangle[p_, mesh_] := Module[
  {tri, \[Xi], ans, distance, minDistance}, 
  ans = {-1, {0, 0}}; 
  minDistance = 10^10; 
  Do[
    tri = mesh["nodes"][[mesh["elems"][[i]]]]; 
    \[Xi] = Chop[IsoparametricCoordinates[p, tri]]; 
    inside = 0 <= \[Xi][[1]] <= 1 && 
             0 <= \[Xi][[2]] <= 1 && 
             0 <= 1 - (\[Xi][[1]] + \[Xi][[2]]) <= 1; 
    If[inside, ans = {i, \[Xi]}; Break[]]; 
    distance = Norm[p - Mean[tri]]; 
    If[distance < minDistance, ans = {i, \[Xi]}; minDistance = distance]
    , 
    {i, 1, mesh["nelems"]}
  ]; 
  ans
]

InterpolationMatrix[coarse_, fine_] := Module[{tri, j, \[Xi], entries},
  entries = Flatten[Table[
    {j, \[Xi]} = FindNearestTriangle[fine["nodes"][[i]], coarse];
    tri = coarse["elems"][[j]];
    {
      {i, tri[[1]]} -> 1 - \[Xi][[1]] - \[Xi][[2]],
      {i, tri[[2]]} -> \[Xi][[1]],
      {i, tri[[3]]} -> \[Xi][[2]]
    },
    {i, 1, fine["nnodes"]}
  ]];
  SparseArray[entries]
]

Note: this FindNearestTriangle[] function is a naive $O(n^2)$$O(n^2)$ 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 ${\mathbf{A}}^{2h}:={\mathbf{R}}^{2h}{\mathbf{A}}^h{\mathbf{I}}^{2h}$$\bold{A}^{2h} \coloneqq \bold{R}^{2h} \bold{A}^{h} \bold{I}^{2h}$

   ➕ Accuracy: this choice has ideal spectral properties

   ➖ Memory Use: requires forming ${\mathbf{A}}^{2h},{\mathbf{A}}^{4h},...$$\bold{A}^{2h}, \bold{A}^{4h}, ...$ 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 ${\mathbf{A}}^{2h},{\mathbf{A}}^{4h},...$$\bold{A}^{2h}, \bold{A}^{4h}, ...$ 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)$$(-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:

(* 
i0 is the depth in the hierarchy, 
k is the number of smoothing iterations at each stage,
  II[[i]] is the interpolation operator at depth i
   A[[i]] is the coefficient matrix at depth i
invD[[i]] is the inverse of the diagonal of A[[i]]
*)
Vcycle[xh0_, rh_, i0_, k_] := Module[{xh = xh0, x2h, r2h},
  If[i0 < nlevels,
    xh = Jacobi[A[[i0]], invD[[i0]], xh, rh, k];
    r2h = Transpose[II[[i0]]] . (rh - A[[i0]] . xh);
    xh += II[[i0]] . Vcycle[r2h*0, r2h, i0 + 1, k];
    xh = Jacobi[A[[i0]], invD[[i0]], xh, rh, k];
    ,
    xh = invA . rh;
  ];
  xh
]

FMGcycle[rh_, i0_, k_] := Module[{xh, x2h, r2h},
  If[i0 < nlevels,
    r2h = Transpose[II[[i0]]] . rh;
    x2h = FMGcycle[r2h, i0 + 1, k];
    xh = II[[i0]] . x2h;
    xh = Vcycle[xh, rh, i0, k]
    ,
    xh = Vcycle[rh*0, rh, i0, k]
  ];
  xh
]

As Fixed-Point iterations

Fixed point iteration was the original iterative way of solving $\mathbf{A}x=b$$\bold{A} x = b$. Starting with some $x_0$$x_0$, subsequent iterates are determined by

for some choice of matrix $\mathbf{M}$$\bold{M}$ (for example, Jacobi method takes $\mathbf{M}=\text{diag}(\mathbf{A}{)}^{-1}$$\bold{M} = \text{diag}(\bold{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 "$\mathbf{M}$$\bold{M}$" corresponding to a 5-deep V cycle with $k$$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 ${\mathbf{A}}^{2h},{\mathbf{A}}^{4h},...$$\bold{A}^{2h}, \bold{A}^{4h}, ...$

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,\mathbf{A}b,{\mathbf{A}}^{2}b,{\mathbf{A}}^{3}b...{\mathbf{A}}^{n}b\}$$\{b, \bold{A}b, \bold{A}^2b, \bold{A}^3 b ... \bold{A}^n b \}$ 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 ($\kappa \approx 30,000,000$$\kappa \approx 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 $\mathbf{M}$$\bold{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.

Conclusions

For this example problem, some key takeaways were that:

• FMG cycles perform better than V cycles, and aren't much more expensive

• Multigrid cycles were more effective as a preconditioner to cg than as a direct fixed point iteration

• Option 2 was less effective than option 1 (but still very good) for V cycles

• Options 1 and 2 were practically equivalent in FMG cycles

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

• How to enforce constraints (e.g. essential boundary conditions, contact) on the coarse systems?

• How well does unstructured multigrid (both options 1 and 2) work in:

  • nearly-incompressible elasticity / saddle point problems

  • problems with strong material anisotropy / heterogeneity

  • problems discretized with Nedelec / Raviart-Thomas elements