In this dissertation, we develop geometric multigrid methods for the finite element approximation of flow problems (e:g:, Stokes, Darcy and Brinkman models) in highly heterogeneous porous media. Our method is based on H^(div)-conforming discontinuous Galerkin methods and the Arnold-Falk-Winther (AFW) type smoothers. The main advantage of using H^(div)-conforming elements is that the discrete velocity field will be globally divergence-free for incompressible fluid flows. Besides, the smoothers used are of overlapping domain decomposition Schwarz type and employ a local Helmholtz decomposition. Our flow solvers are the combination of multigrid preconditioners and classical iterative solvers. The proposed preconditioners act on the combined velocity and pressure space and thus does not need a Schur complement approximation. There are two main ingredients of our preconditioner: first, the AFW type smoothers can capture a meaningful basis on local divergence free subspace associated with each overlapping patch; second, the grid operator does not increase the divergence from the coarse divergence free subspace to the fine one as the divergence free spaces are nested. We present the convergence analysis for the Stokes' equations and Brinkman's equations ( with constant permeability field ), as well as extensive numerical experiments. Some of the numerical experiments are given to support the theoretical results. Even though we do not have analysis work for the highly heterogeneous and highly porous media cases, numerical evidence exhibits strong robustness, efficiency and unification of our algorithm.