We will introduce a first-order system least squares (FOSLS) formulation of the diffusion equation with discontinuous coefficient. We discretize this equation using a standard finite element space augmented by singular basis functions. These singular basis functions enable us to resolve the flux without deterioration of discretization order in the neighborhood of singularities. Such singularities of the flux are typically induced by the geometry of the discontinuity of the porosity coefficient, or by reentrant corners. We introduce a multilevel solver for the discrete problem that utilizes singular basis functions on coarser levels, without compromising the efficiency of the solver. This research is supported by hte Department of Energy, under contract W-7405-ENG-36, LA-UR-00-575.