A computer-implemented method for performing simulation of a reservoir represented by a model of a reservoir pore system having a plurality of matrix nodes (203, 302, 704a-c), a plurality of fracture nodes (201, 304, 702a-f), at least one fracture zone and at least one fracture-free zone , and wherein the matrix nodes represent pore systems comprising high pore volume and low conductivity and wherein fracture nodes represent pore systems comprising low pore volume and high conductivity, the method comprising: for fracture subgrids (202a-202d) associated with matrix subgrids (200a-200d), characterizing each of the plurality of matrix nodes as dual matrix nodes or non-dual matrix nodes, wherein the fracture subgrid comprises a subset of fracture nodes and the matrix subgrid comprises a subset of matrix nodes, wherein associated subgrids represent subgrids with the same physical location but different porosity type, wherein dual matrix nodes are matrix nodes connected to fracture nodes, and wherein non-dual matrix nodes are matrix nodes connected to other matrix nodes, characterised by : introducing virtual fracture nodes in areas of the fracture subgrids with missing fracture nodes, wherein virtual fracture nodes are dual matrix nodes without a fracture node associated therewith, wherein the virtual fracture nodes are introduced through matrix-to-matrix connections or fracture-to-matrix connections; eliminating dual matrix nodes or dual matrix blocks and updating a coefficient matrix and Right Hand Side, RHS, vector for the model accordingly, wherein dual matrix blocks are comprised of one or more subdivided dual matrix nodes, wherein the subdivided dual matrix nodes are smaller in bulk volume compared to fracture nodes and matrix nodes; determining in-fill patterns between fracture nodes and virtual fracture nodes with matrix-to-matrix connections, the in-fill patterns resulting from the elimination; where virtual fracture nodes are introduced, merging the virtual fracture nodes to the fracture subgrid and updating the coefficient matrix and RHS vector for the model accordingly; solving the resultant linear system using a linear solver ; and backsolving for the eliminated dual matrix nodes, thereby simulating the reservoir whereby, based upon the reservoir simulation, a drilling plan may be implemented. introducing virtual fracture nodes in areas of the fracture subgrids with missing fracture nodes, wherein virtual fracture nodes are dual matrix nodes without a fracture node associated therewith, wherein the virtual fracture nodes are introduced through matrix-to-matrix connections or fracture-to-matrix connections; eliminating dual matrix nodes or dual matrix blocks and updating a coefficient matrix and Right Hand Side, RHS, vector for the model accordingly, wherein dual matrix blocks are comprised of one or more subdivided dual matrix nodes, wherein the subdivided dual matrix nodes are smaller in bulk volume compared to fracture nodes and matrix nodes; determining in-fill patterns between fracture nodes and virtual fracture nodes with matrix-to-matrix connections, the in-fill patterns resulting from the elimination; where virtual fracture nodes are introduced, merging the virtual fracture nodes to the fracture subgrid and updating the coefficient matrix and RHS vector for the model accordingly; solving the resultant linear system using a linear solver ; and backsolving for the eliminated dual matrix nodes, thereby simulating the reservoir whereby, based upon the reservoir simulation, a drilling plan may be implemented. <EOS>