In computer implemented simulation of fluid characteristics of a subsurface reservoir with a reservoir simulator, the fluid characteristics being governed by fluid flow, material balance and pressure-volume-temperature relationships, the reservoir being partitioned into a number of subdomains, each subdomain representing a finite volume of the reservoir containing multiphase, multicomponent fluids and being organized into a group of grid cells having data regarding the fluids, the reservoir simulation being performed in a cluster computer formed of at least one master node and a plurality of processor nodes, with each processor node being assigned a subdomain of the partitioned reservoir, a computer implemented method to simulate by Newtonian iterations fluid characteristics parameters of interest for the grid cells, a computer implemented method comprising the steps of: arranging the data for each of the reservoir subdomains into a Jacobian matrix A composed of: a matrix P composed of the block diagonal of the subdomain and non-zero data blocks adjacent the block diagonal, a matrix E composed of non-zero data blocks not included in the block diagonal matrix P, the matrix E being composed of a matrix E' of interior grid blocks within the subdomain and a matrix E" of boundary grid cell blocks with grid cells of adjacent subdomains in the reservoir; performing parallel approximate solutions of the fluid characteristics for the grid blocks of each of the subdomains in the assigned processor nodes using a chosen local preconditioner for the subdomain; updating full system residuals for the Jacobian matrix for the reservoir based on the parallel approximate solutions; restricting the full system residuals into a global reduced space according to a global reduced space matrix; performing a parallel approximate solution of the full system residuals using a global reduced space preconditioner; forming a combined approximate solution update in the assigned processor nodes of the parallel approximate solutions of each of the subdomains and the parallel approximate solution of the full system residuals to form a solution vector; updating the full system residuals using the solution vector; applying an iterative solver to the full system to obtain measures of the fluid material balances of each composition for the grid cells of the reservoir; testing the measures of the fluid material balances of each composition obtained by applying the iterative solver to determine if convergence is achieved; and if convergence is not achieved, returning to the step of performing parallel approximate solutions; or if convergence is achieved, returning to the reservoir simulator. arranging the data for each of the reservoir subdomains into a Jacobian matrix A composed of: a matrix P composed of the block diagonal of the subdomain and non-zero data blocks adjacent the block diagonal, a matrix E composed of non-zero data blocks not included in the block diagonal matrix P, the matrix E being composed of a matrix E' of interior grid blocks within the subdomain and a matrix E" of boundary grid cell blocks with grid cells of adjacent subdomains in the reservoir; performing parallel approximate solutions of the fluid characteristics for the grid blocks of each of the subdomains in the assigned processor nodes using a chosen local preconditioner for the subdomain; updating full system residuals for the Jacobian matrix for the reservoir based on the parallel approximate solutions; restricting the full system residuals into a global reduced space according to a global reduced space matrix; performing a parallel approximate solution of the full system residuals using a global reduced space preconditioner; forming a combined approximate solution update in the assigned processor nodes of the parallel approximate solutions of each of the subdomains and the parallel approximate solution of the full system residuals to form a solution vector; updating the full system residuals using the solution vector; applying an iterative solver to the full system to obtain measures of the fluid material balances of each composition for the grid cells of the reservoir; testing the measures of the fluid material balances of each composition obtained by applying the iterative solver to determine if convergence is achieved; and if convergence is not achieved, returning to the step of performing parallel approximate solutions; or if convergence is achieved, returning to the reservoir simulator.