It is critical to develop numerical simulators that can robustly and efficiently predict production from unconventional reservoirs. In this paper, we describe in detail the formulation, flexibility and simulation results of our recently developed shale gas simulator. The first component of this work includes Discrete Fracture Modeling (DFM) of hydraulic fractures. To capture the flow transients near fracture stages and horizontal well, we develop a gridding algorithm based on Centroidal Voronoi Tessellation (CVT). To represent the hydraulic fractures, we populate nodes around fracture stages based on the refinement that is required by the user. The populated nodes are viewed as constraints and are kept stationary during the subsequent iterations. For matrix region we use a two-step algorithm to generate the mesh. In the first step, we tessellate the domain by minimizing the forces along the bars between adjacent nodes, analogues to mass-spring or truss systems. In the next step, we perform an iterative gird optimization to generate Centroidal Voronoi Tessellation. The second component of the paper targets representation of a wide range of physical phenomena including gas slippage, gas adsorption/desorption, and high velocity non-Darcy flow of gas. Our formulation is fully implicit and allows for flow of up to two phases (water and gas). To avoid manual calculation of the Jacobian corresponding to the coupled nonlinear system of equations, we utilize Automatic Differentiation (AD), allowing us to readily incorporate the nonlinear governing equations at nearly no extra cost. We perform sensitivity studies to investigate the effect of different parameters on production from shale gas plays as well as on computational efficiency of the simulator. (C) 2015 Elsevier Ltd. All rights reserved.