The displacement formulation of linear elasticity yields a system of equations in which the grad-div operator is the dominant term in the incompressible limit. A similar operator also arises when using least-squares techniques to pose the linear Boltzmann equation as a minimization problem. When using a standard FE discretization (using linear FEs) and a standard multigrid method (using pointwise smoothing) for either system, non-robust convergence results will be observed in the multigrid method as the grad-div operator becomes dominant. The prescription to overcoming this insufficient performance is to use a FE space with a sizable divergence-free subspace and the proper relaxation operator. In this talk, we investigate using different orders of polynomials to discretize the problem and a group smoothing approach to reduce divergence-free error in the multigrid routine.