Parallel Domain Decomposition Solvers for Large-Scale Reservoir Simulation

Contacts: Mary F. Wheeler [email], Marcelo Ramé [email], Carol Woodward [email].

Domain decomposition methods are a major area of contemporary researchin numerical analysis of partial differential equations. They providerobust, parallel and scalable preconditioned iterative methods for thelarge-scale linear systems arising from discretization of continuousproblems by finite elements, finite differences or spectral methods.

In reservoir simulation, large systems of linear or nonlinearequations have to be solved at each time step, with one or more unknownsper grid block, depending on the formulation of the discrete problem. These linear systems or the linearization of nonlinear systems resultingfrom Newton’s method give rise to sparse, banded matrices amenable tosolution by iterative methods. The spatial variability of rock propertiesin a reservoir, as well as flow conditions can yield mobility coefficientsjumping four to five orders of magnitude over short distances. As aresult of this, the associated algebraic system’s matrices are highlyill-conditioned, thus requiring the use of robust preconditioners fortheir iterative solution.

Two such preconditioners for the conjugate gradient method, which showthe required robustness as well as excellent parallel scaling propertieswere implemented for the three-dimensional, parallel simulatorRice-UTCHEM. The capabilities of this code are presented elsewhere inthese pages. The numerical discretization (IMPES) requires the solutionof a linear system of pressure unknowns at each time step. All othervariables of the model are updated explicitly. The block-centeredfinite-difference discretization (equivalent to the lowest order mixedfinite elements) gives rise to a pressure linear system which is symmetricpositive definite.

Both the Balancing Domain Decomposition (BDD) and the OverlappingAdditive Schwarz (OAS) preconditioners for the conjugate gradient methodare highlighted on this page. The BDD preconditioner is a nonoverlappingmethod that solves for the interface unknowns between subdomains asprimary unknowns (details of the method are given in [1]). The OASpreconditioner is an overlapping preconditioner and is implemented herewith minimal overlap, i.e., thickness = 1 grid-block (details are given in[2]). Both domain decomposition preconditioners are compared against thediagonal or Jacobi preconditioner in the figure shown. The matrix andvector of right hand sides were generated by the parallel simulatorRiceUTCHEM for a 3-D case with a mobility contrast of 3 orders ofmagnitude. This is a problem of only modest complexity, both from thestandpoint of the number of unknowns and the distribution of mobilitiesthroughout the reservoir.



The above figure shows, on a log-log scale, the elapsed times required tosolve this problem, on the Touchstone Delta (Intel massively-parallelprototype machine at CalTech), using JCG, BDD-CG and OAS-CG. Decompositions from 8 to 128 subdomains are shown on the graph. The solidblack line, with a negative unit slope, indicates the ideal performance,i.e., doubling the number of processors reduces the time by a half. BothDD methods solve this problem much faster than the Jacobi preconditioner. Additionally, the poor scaling of the JCG is evidenced by the slope of thecorresponding curve, which is under unity in absolute value for all casestried. The DD methods show much better scaling features up to the point(around 100 processors for this problem) where the subdomains become tosmall and interprocessor communication starts to dominate the elapsedtimes.