Multigrid methods have become the optimal solvers for Lattice QCD, first for the Wilson-clover discrete representation and more recently for the domain wall formulation. However, at present, the ensembles with the largest lattices use a third discretization, staggered. At the physical pion mass, propagator inversions with the staggered operator require $\mathcal{O}(10,000)$ iterations using the Conjugate Gradient algorithm. The staggered discretization raises new challenges for a multigrid implementation. It is a first order discretization which in four-dimensions corresponds to four copies of a Dirac fermion in the continuum. Here we report on our investigation into a new geometric adaptive multigrid for staggered fermions both for the normal equations and for a first order projection based on a novel blocking scheme stabilized by a second order gauged Laplacian. Current numerical tests, applied to the two-dimensional staggered representation of the two-flavor Schwinger model, are promising. Further improvements are under way as well as generalizations and tests for four-dimensional QCD using multigrid preconditioned staggered solvers in the QUDA library for multi-GPU architectures.