The objective of this work is to develop a robust and efficient multilevel parallel algorithm to solve regelar composite grids resulting from local sub-gridding of a multi-million-cell base grid. The strategy is to formulate an iterative scheme on a grid-by-grid basis such that the regelar structure of the matrix at each stage is preserved. The multilevel method applied here utilizes an algebraic fill-in method similar to the one proposed in BEPS but has been generalized to the multi-component multi-phase flow context.