Unstructured Multigrid for Elasticity

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.

Interpolation and Restriction

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 ${\mathbf{I}}^{2h},{\mathbf{R}}^{2h}$$\bold{I}^{2h}, \bold{R}^{2h}$ can still work for transferring vector-valued fields between meshes, provided the nodal values are arranged as a $n\times 2$$n\times2$ matrix. However, the "option 1" approach of defining the coarse coefficient matrices by

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 $[u_{1x},u_{1y},u_{2x},u_{2y},...]$$[u_{1x}, u_{1y}, u_{2x}, u_{2y}, ... ]$), define

where ${\mathbf{1}}_d$$\bold{1}_d$ is the identity matrix with $d$$d$ rows and columns and $\otimes$$\otimes$ is the Kronecker product. This lets us update our expression for the coarse elasticity operators to

Numerical Example

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

where $\mathbf{K}$$\bold{K}$ is the stiffness matrix for a plane strain elasticity problem (with $E=100,\nu =0.25$$E=100, \;\nu = 0.25$, linear isotropic) and $\mathbf{M}$$\bold{M}$ is the consistent mass matrix (with $\rho =1$$\rho = 1$). If we form these coefficient matrices from the wrench mesh hierarchy and run the previous multigrid test suite, we get the following

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 $\omega$$\omega$ ("regular" Jacobi is $\omega =1$$\omega = 1$). Choosing $\omega \approx \frac{2}{{\lambda }_{max}({\mathbf{D}}^{-1}\mathbf{A})}$$\omega \approx \frac{2}{\lambda_{max}(\bold{D}^{-1} \bold{A})}$ makes all the eigenvalues of the iteration matrix fall inside the unit disk again. In this next plot, we show the convergence of the different multigrid options for two choices of $\omega$$\omega$: the solid lines represent the conservative choice of $\omega =0.5$$\omega = 0.5$ at all stages, while the dashed lines use the optimal (but more expensive to determine) choice ${\omega }^h=\frac{2}{{\lambda }_{max}(\mathbf{D}{{}^h}^{-1}{\mathbf{A}}^h)}$$\omega^{h} = \frac{2}{\lambda_{max}(\bold{D}{^h}^{-1} \bold{A}^{h})}$.

So, both choices of $\omega$$\omega$ make multigrid-as-a-preconditioner effective again. If we know the maximum eigenvalue of ${\mathbf{D}}^{-1}\mathbf{A}$$\bold{D}^{-1} \bold{A}$, then the optimal choice of $\omega$$\omega$ does help convergence, but only by 10-15%. For a problem with multiple linear solves that use the same coefficient matrix, it may be worthwhile to estimate that maximum eigenvalue up front, in order to accelerate convergence of subsequent solves.

Comparison With Hypre

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.

Conclusion

After making two small modifications (i.e. $\hat{\mathbf{I}}$$\hat{\bold{I}}$ vs $\mathbf{I}$$\bold{I}$, and moving to $\omega -$$\omega-$weighted Jacobi smoothing), multigrid-as-a-preconditioner works well for this plane strain elasticity problem. Much like the Poisson problem we find that:

• more smoothing per stage $\Rightarrow$$\Rightarrow$ faster convergence per iteration, but more expensive iterations

• 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 ${\lambda }_{max}({\mathbf{D}}^{-1}\mathbf{A})$$\lambda_{max}(\bold{D}^{-1} \bold{A})$ to make the weighted smoothing robust. Another alternative is using a different smoother whose iteration matrix eigenvalues naturally lie in the unit disk.