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
, where is some quantity, and represents the element size. The analogous quantities on the coarser and coarser meshes will be referred to as 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.
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
The interpolation operators can be defined unambigously for continuous finite element fields by recognizing that the mesh nodes themselves should satisfy
where
The isoparametric coordinates are used to evaluate that element's shape functions, which become the values of the nonzero entries in that row of
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
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 naiveimplementation. For large meshes, you would want to use something like a BVH to accelerate this sort of spatial query.
The last thing we need is a way to represent the physics on the coarse meshes. There are two natural ways to do this:
Use the interpolation/restriction operators to form
➕ Accuracy: this choice has ideal spectral properties
➖ Memory Use: requires forming
➖ Complexity: sparse matrix assembly is expensive and complicated
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
➖ 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
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).
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
]
Fixed point iteration was the original iterative way of solving
for some choice of matrix
we see that the residual drops quickly for the first few iterates, but convergence slows down and stagnates. Compare that to "
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
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.
Krylov methods are the more modern way of solving systems of equations iteratively. They use the subspace spanned by cg
) to our system of equations (
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
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.
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:
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