The initial writeup looked at applying multigrid to a Poisson problem defined on an unstructured mesh. This document will apply the same ideas to plane strain elasticity.

The interpolation and restriction operators change slightly with elasticity problems, since each node has more than one degree of freedom. So, the dimensions of the coefficient matrix are bigger by a factor of 2 (since this is a 2D problem). The original definitions of

doesn't quite work since there's a size mismatch between the bigger coefficient matrix and the interpolation/restriction operators. For this case (assuming the displacement vector is ordered

where

For these examples, we'll use a coefficient matrix with the following form

where

Most of the approaches that worked for the Poisson problem are either not converging, or are converging to the wrong answer. The V cycles are making the residual smaller, but progress is incredibly slow.

What is going on?

Regular Jacobi method isn't working for this elasticity problem out of the box, and all of our multigrid cycles are built on Jacobi smoothing. If we look at the largest eigenvalues of the Jacobi iteration matrix from the finest mesh, we see

The presence of eigenvalues greater than 1 means that those corresponding mode shapes get amplified with each Jacobi iteration. This can be dealt with by modifying the smoothing to be a "weighted Jacobi iteration" instead

Here, we have to choose

So, both choices of

This same problem was set up with mfem + hypre, using the BoomerAMG preconditioner with default settings. As I understand it, there are two ways in mfem to configure BoomerAMG to use its elasticity-specific variant:

by calling

`mfem::HypreBoomerAMG::SetSystemsOptions()`

by calling

`mfem::HypreBoomerAMG::SetElasticityOptions()`

It's unclear which of these is preferred, so I include both in the following comparison. The console output from hypre reveals that mfem's default BoomerAMG configuration seems to be a single V cycle, with 1 smoothing iteration per level, so that's what we'll compare against.

In this example, BoomerAMG's preconditioners both perform poorly, failing to make progress until ~600 iterations in. Our implementation of multigrid-as-a-preconditioner converges monotonically and takes a fraction of the iterations that BoomerAMG does.

The nonconvergence of BoomerAMG for a simple elasticity problem is concerning, but this issue has been brought to the attention of mfem and hypre developers several times without resolution. So, it appears it might not be just user error. For anyone interested in debugging this problem, logs of hypre's console output for the two preconditioners are available here:

hypre_output_set_systems_options

hypre_output_set_elasticity_options

Reading these log files, it sounds like hypre may be performing smoothing operations at the coarsest level, rather than factorizing the matrix and solving. If this is true, it seems like a critical oversight, as the cost of factorizing the coarsest matrix is totally insignificant and it makes a huge difference in the effectiveness of multigrid methods as a whole. For example, if I replace the linear solve at the coarsest level by smoothing operations, the multigrid methods take an order of magnitude more iterations to arrive at the solution:

I was unable to measure any significant difference in time per iteration between the two approaches (i.e. linear solve vs. smoothing at the coarsest level), so the only metric that matters in this comparison is iteration count.

After making two small modifications (i.e.

more smoothing per stage

faster convergence per iteration, but more expensive iterations$\Rightarrow $ FMG converges in fewer iterations than V-cycles

both interpretations of coarse coefficient matrices are effective

From an implementation standpoint, it would be nice to find a procedure for either bounding, estimating, or quickly computing