A computer-implemented method of performing a reservoir simulation to simulate a property of at least one fluid in a fluid-containing system such as a reservoir, comprising: generating a data representation of a computational mesh (200, 400) representing a reservoir in a storage system, wherein the data representation comprises an interconnection weight that represents the magnitude of an interconnection between each of a plurality of computational cells (202, 402) in the computational mesh, and wherein: the computational mesh is partitioned into subdomains (I, II, III, IV) and each subdomain is assigned to a different computing unit or each subdomain is assigned to multiple computing units or multiple subdomains are assigned to a single computing unit, the data representation includes a Jacobian matrix in which each diagonal element represents one of the plurality of computational cells in the computational mesh, and the interconnection weights are generated between each of the plurality of computational cells from off-diagonal elements in the Jacobian matrix; determining a threshold value, wherein determining the threshold value comprises: calculating a maximum interconnection weight for interconnections of each of the plurality of computational cells to all adjoining computational cells, identifying the largest value of the maximum interconnection weight for all of the plurality of computational cells, identifying the smallest value of the maximum interconnection weight for all of the plurality of computational cells, and selecting the threshold value to lie between the largest value and the smallest value; comparing the threshold value to each interconnection weight using the formula 1 + d ij + d ji d ij d ji a ij ≤ σ, wherein d ij and d ji represent the diagonal elements of the computational cells i and j in the Jacobian matrix, σ represents the threshold value, and a ij represents transmissibility; setting any interconnection weight that is equal to or less than the threshold value to zero; and eliminating any rows in the Jacobian matrix that have no remaining off-diagonal elements to generate a lower level matrix. the computational mesh is partitioned into subdomains (I, II, III, IV) and each subdomain is assigned to a different computing unit or each subdomain is assigned to multiple computing units or multiple subdomains are assigned to a single computing unit, the data representation includes a Jacobian matrix in which each diagonal element represents one of the plurality of computational cells in the computational mesh, and the interconnection weights are generated between each of the plurality of computational cells from off-diagonal elements in the Jacobian matrix; determining a threshold value, wherein determining the threshold value comprises: calculating a maximum interconnection weight for interconnections of each of the plurality of computational cells to all adjoining computational cells, identifying the largest value of the maximum interconnection weight for all of the plurality of computational cells, identifying the smallest value of the maximum interconnection weight for all of the plurality of computational cells, and selecting the threshold value to lie between the largest value and the smallest value; comparing the threshold value to each interconnection weight using the formula 1 + d ij + d ji d ij d ji a ij ≤ σ, wherein d ij and d ji represent the diagonal elements of the computational cells i and j in the Jacobian matrix, σ represents the threshold value, and a ij represents transmissibility; setting any interconnection weight that is equal to or less than the threshold value to zero; and eliminating any rows in the Jacobian matrix that have no remaining off-diagonal elements to generate a lower level matrix. calculating a maximum interconnection weight for interconnections of each of the plurality of computational cells to all adjoining computational cells, identifying the largest value of the maximum interconnection weight for all of the plurality of computational cells, identifying the smallest value of the maximum interconnection weight for all of the plurality of computational cells, and selecting the threshold value to lie between the largest value and the smallest value; comparing the threshold value to each interconnection weight using the formula 1 + d ij + d ji d ij d ji a ij ≤ σ, wherein d ij and d ji represent the diagonal elements of the computational cells i and j in the Jacobian matrix, σ represents the threshold value, and a ij represents transmissibility; setting any interconnection weight that is equal to or less than the threshold value to zero; and eliminating any rows in the Jacobian matrix that have no remaining off-diagonal elements to generate a lower level matrix.