Subsurface rocks are naturally porous and may contain numerous fractures that intersect in complex configurations. Fractures are, in general, permeable pathways and are responsible for most of the flow. They are classically modeled as structures of co-dimension one in the surrounding rock. Due to the porosity of the rock, a 3D flow also occurs in the rock, and this flow is coupled to the 2D flow in the fractures. The aim of this work is to efficiently simulate single-phase flow in large-scale fractured porous media. Meshing this geometry is a challenge that we have taken up with the development of a dedicated surface mesher, MODFRAC, and the use of the efficient volume mesher, GHS3D. Due to the severe constraints induced by the fractures, the mesh may contain low-quality elements. The mixed hybrid finite element (MHFEM) is convenient for discretizing this coupled flow problem, as the face unknowns ease the coupling between the 2D (in the fractures) and 3D (in the rock) flows. In addition, the MHFEM leads to a symmetric positive definite linear system with traces of hydraulic head as unknowns, the velocity then being calculated locally, in a post-processing step. However, the linear system may contain millions of unknowns. Direct solvers are no longer an option due to excessive RAM consumption. Moreover, the matrix of the linear system is very poorly conditioned due to the permeability contrasts between rock and fractures, as well as the low-quality elements present in the mesh. In this presentation, we will show how to set up the data to use the domain decomposition HPDDM library, implemented in PETSc, as preconditioner. We will present several examples demonstrating the excellent performance obtained with HPDDM on highly heterogeneous and large-scale fractured porous media, containing up to 700k fractures (242M of unknowns).