Parallel Implementation of a Two-level Algebraic ILU(k)-based Domain Decomposition Preconditioner ^†

ABSTRACT

We discuss the parallel implementation of a two-level algebraic ILU(k)-based domain decomposition preconditioner using the PETSc library. We present strategies to improve performance and minimize communication among processes during setup and application phases. We compare our implementation with an off-the-shelf preconditioner in PETSc for solving linear systems arising in reservoir simulation problems, and show that for some cases our implementation performs better.

Keywords:
Two-level preconditioner; domain decomposition; Krylov methods; linear systems; parallelism; PETSc

RESUMO

Discutimos a implementação paralela de um precondicionador algébrico de decomposição o de domínios em dois níveis baseado em ILU(k), utilizando a biblioteca PETSc. Apresentamos estratégias para melhorar a performance, minimizando a comunicação entre processos durante as fases de construção e de aplicação. Comparamos, na solução de sistemas lineares provenientes de problemas de simulação de reservatórios, a nossa implementação com um precondicionador padrão do PETSc. Mostramos que, para alguns casos, nossa implementação apresenta um desempenho superior.

Palavras-chave:
Precondicionador de dois níveis; decomposição de domínio; métodos de Krylov; sistemas lineares; paralelismo; PETSc

1 INTRODUCTION

This paper discusses the formulation and the parallel implementation of an algebraic ILU(k)-based two-level domain decomposition preconditioner first introduced in².

In this work we present and discuss details of the implementation using the MPI-based PETSc suite³, a set of data structures and routines for the parallel solution of scientific applications modeled by partial differential equations. We also present results of computational experiments involving matrices from oil reservoir simulation. We have tested with different number of processes and compared the results with the default PETSc preconditioner, block Jacobi, which is a usual option in the oil industry.

The multilevel preconditioner has been an active research area for the last 30 years. One of the main representatives of this class is the algebraic multigrid (AMG) method^{19), (22), (30), (31}, when used as a preconditioner rather than as a solver. It is widely used in oil reservoir simulation as part of the CPR (constrained pressure residual) preconditioner^6),(33. Despite its general acceptance there is room for new alternatives, as there are problems where AMG can perform poorly, see, for instance¹³. Among these alternatives we find the two-level preconditioners. There are many variations within this family of preconditioners, but basically we can discern at least two subfamilies: algebraic^{1), (4), (12), (14), (18), (20), (23), (28), (32} and operator dependent^{15), (16), (17}. Within the operator dependent preconditioners we should highlight the spectral methods^{21), (25}.

Incomplete LU factorization ILU(k)²⁶ has long been used as a preconditioner in reservoir simulation (an ingenious parallel implementation is discussed in¹¹. Due to the difficulty in parallelizing ILU(k), it is quite natural to combine ILU(k) and block-Jacobi, so much so that this combination constitutes PETSc’s default parallel preconditioner³. The algorithm proposed in², whose parallel implementation we discuss in this article, seeks to combine the use of (sequential) ILU(K) with two ideas borrowed from domain decomposition methods: (i) the introduction of an interface that connects subdomains, allowing, as opposed to block-Jacobi, for the interaction between subdomains to be taken into account, and (ii) the introduction of a second level, associated to a coarse version of the problem, that speeds up the resolution of low frequency modes. These improvements come at the cost of greater communication, requiring a more involved parallel implementation.

The main contribution of the two-level preconditioner proposed in² is the fine preconditioner, as the coarse level component is a quite simple one. Accordingly, the main contribution in our article is the development of a careful parallel implementation of that fine part; nonetheless, we also take care of the coarse part by proposing low-cost PETSc-based parallel codes for its construction and application. Besides the parallel implementation, we present and discuss a set of performance tests of this preconditioner for solving synthetic and real-world oil reservoir simulation problems.

The paper is organized as follows. Section 2 introduces the notation used throughout the paper. Section 3 describes in detail the proposed two-level algebraic ILU(k)-based domain decomposition preconditioner, which we call iSchur. Its parallel implementation is detailed in Section 4, including the communication layout and strategies to minimize data transfer between processes, while results of performance experiments and comparisons with another preconditioner are presented and discussed in Section 5. The conclusion is drawn in Section 6 along with some future work directions.

2 NOTATION

In this section, we introduce some notation that will be necessary to define iSchur.

Consider the linear system of algebraic equations

$Ax=b$ (2.1)

arising from the discretization of a system of partial differential equations (PDEs) modeling the multiphase flow in porous media by a finite difference scheme with n gridcells. We denote by n _dof the number of degrees of freedom (DOFs) per gridcell, so that A is a matrix of dimension n _dof n × n _dof n. For instance, for the standard black-oil formulation, see²⁴, the DOFs are oil pressure, oil saturation and water saturation, so that n _dof = 3.

It will be convenient to think of A as a block-matrix:

$A=\left[\begin{array}{cccc}{a}_{11}& a_{12}& \cdots & a_{1n}\\ a_{21}& a_{22}& \cdots & a_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ a_{n1}& a_{n2}& \cdots & a_{nn}\end{array}\right],$

where a block, denoted a _ij , is a matrix of dimension n _dof × n _dof.

The block-sparsity pattern of A is determined by the stencil of the finite-difference scheme used for the discretization of the PDEs. Figure 1 depicts a bidimensional domain discretized by a 3×3 mesh and the sparsity pattern of the associated matrix for the black-oil model, assuming a five-point finite-difference stencil. We say that two gridcells i and j are neighbors if either one of the blocks a _ij or a _ji is not null. For the five-point stencil, gridcells are neighbors when they share an edge.

Figure 1:
Discretization grid and associated matrix, assuming three unknowns per gridcell and a five-point stencil.

The gridcells and their respective indices are identified, in a way that Ω denotes either the domain of the PDE or the set of indices {1,2,…,n}. We introduce a disjoint partition of Ω, i.e., we break Ω ={1,2,…,n} into P pieces Ω₁, Ω₂,..., Ω_p , in such a way that each gridcell belongs to exactly one subdomain Ω_k , see Figure 2. More precisely, we define

${\left\{{\Omega }_J\right\}}_{1\le J\le P}\begin{array}{c}\end{array}\text{s}\text{.t}\text{.}\begin{array}{c}\end{array}{\displaystyle \underset{J=1}{\overset{P}{\cup }}{\Omega }_J=\Omega \begin{array}{c}\end{array}\text{and}\begin{array}{c}\end{array}{\Omega }_I\cap {\Omega }_J=}\cancel{0}\forall I\ne J.$ (2.2)

Figure 2:
2D Domain partitioned into 4 subdomains.

We note that there are gridcells that, while belonging to one subdomain, have neighbors in other subdomains. Our domain decomposition approach requires the definition of disconnected subdomains, that can be dealt with in parallel. For that sake, we define a separator set, called interface and denoted by Γ, and (disconnected) subdomain interiors ${\Omega }_J^{Int}$. A gridcell j is said to belong to Γ if it is in a subdomain Ω_J while being neighbor of at least one gridcell in another subdomain with greater index, i.e., Ω_K with K > J:

$\Gamma =\left\{j\in \Omega |\exists k,J,K\begin{array}{c}\end{array}\text{s}\text{.t}\text{.}\begin{array}{c}\end{array}j\in {\Omega }_J,k\in {\Omega }_K,K>J,\left(a_{jk}\ne 0\begin{array}{c}\end{array}or\begin{array}{c}\end{array}{a}_{kj}\ne 0\right)\right\}.$ (2.3)

We now define the subdomain interior ${\Omega }_J^{Int}$ as the portion of Ω_J not in Γ:

${\Omega }_J^{Int}={\Omega }_J-\Gamma ,$ (2.4)

see Figure 2. These definitions are just enough to ensure that if gridcell j is in ${\Omega }_J^{Int}$ and gridcell k is in ${\Omega }_K^{Int}$, with J ≠ K, then they are not neighbors.¹ Indeed, to fix the notation, assume J < K. Since $j\in {\Omega }_J^{Int}\subset {\Omega }_J$ and $k\in {\Omega }_K^{Int}\subset {\Omega }_K$, if j and k were neighbors, j would be in Γ by definition (2.3) and therefore not in ${\Omega }_J^{Int}$.

We now define the local interface Γ_J associated with each subdomain Ω_J as the intersection of Γ and Ω_J , or equivalently,

${\Gamma }_J=\left\{j\in {\Omega }_J|\left(\exists K>J\begin{array}{c}\end{array}and\begin{array}{c}\end{array}\exists k\in {\Omega }_K\right)\text{s}\text{.t}\text{.}\left(a_{jk}\ne 0\begin{array}{c}\end{array}or\begin{array}{c}\end{array}{a}_{kj}\ne 0\right)\right\}.$ (2.5)

Notice that {Γ_J }_{1≤ J ≤ p} form a disjoint partition of Γ. See Figure 2.

Finally, we define extended subdomains ${\overline{\Omega }}_J$ and extended local interfaces ${\overline{\Gamma }}_J$, which incorporate the portions of the interface connected to the subdomain interior ${\Omega }_J^{Int}$. We define ${\overline{\Omega }}_J$ as

${\overline{\Omega }}_J={\Omega }_J\cup \left\{k\in \Gamma |\exists j\in {\Omega }_J^{Int}\begin{array}{c}\end{array}\text{s}\text{.t}\text{.}\begin{array}{c}\end{array}\left(a_{jk}\ne 0\begin{array}{c}\end{array}or\begin{array}{c}\end{array}{a}_{kj}\ne 0\right)\right\}.$ (2.6)

and ${\overline{\Gamma }}_J$ as its restriction to Γ, i.e., ${\overline{\Gamma }}_J=\Gamma \cap {\overline{\Omega }}_J$, see Figure 2.

Notice that ${\Gamma }_J\subset {\overline{\Gamma }}_J\subset \Gamma$. We point out that ${\overline{\Gamma }}_J$ is the result of augmenting Γ_J with the gridcells of the interface Γ that are neighbors of gridcells in Ω_j . We refer to ${\overline{\Gamma }}_J$ as an extended interface, see Figure 2.

If the equations/variables are reordered, starting with the ones corresponding to ${\Omega }_1^{Int}$, followed by the other ${\Omega }_J^{Int}$ and finally by Γ, then A has the following block-structure:

$A=\left[\begin{array}{cccc}{A}_{11}& & & A_{1\Gamma }\\ & \ddots & & ⋮\\ & & A_{PP}& A_{P\Gamma }\\ A_{\Gamma 1}& \cdots & A_{\Gamma P}& A_{\Gamma \Gamma }\\ & & & \end{array}\right],$ (2.7)

where the submatrices A _JJ contain the rows and columns of A associated with ${\Omega }_J^{Int}$, A _ΓΓ the ones associated with Γ, A _JΓ the rows associated with ${\Omega }_J^{Int}$ and the columns associated with Γ, and A _ΓJ the rows associated with Γ and the columns associated with ${\Omega }_J^{Int}$. Therefore, denoting by |S| the number of elements in a set S, the dimensions of A _JJ , A _ΓΓ, A _JΓ and A _ΓJ are, respectively, n _dof n _J × n _dof n _J, n _dof m× n _dof m, n _dof n _J × n _dof m, and n _dof m× n _dof n _J , where $n_J=\left|{\Omega }_J^{Int}\right|$ and m=|Γ|. It is important to notice that, since a _jk = 0 for any $j\in {\Omega }_J^{Int}$ and $k\in {\Omega }_K^{Int}$ with J ≠ K, the submatrices A _JK of rows associated with ${\Omega }_J^{Int}$ and columns associated with ${\Omega }_K^{Int}$ are all null (and therefore omitted in (2.7)). This block-diagonal structure of the leading portion of the matrix (which encompasses most of the variables/equations) allows for efficient parallelization of many tasks, and was the ultimate motivation of all the definitions above. We point out that although not necessary by the method, the implementation discussed in this work assumes that the matrix A is structurally symmetric², in this case the matrices A _JJ and A _ΓΓ are also structurally symmetric. Furthermore, A _JΓ and $A_{\Gamma J}^T$ have the same nonzero pattern.

3 DESCRIPTION OF THE TWO-LEVEL PRECONDITIONER

The goal of a preconditioner is to replace the linear system (2.1) by one of the equivalent ones:

$\begin{array}{l}\left(MA\right)x=Mb\begin{array}{c}\end{array}\left(\text{left preconditioned}\right)\text{or}\\ \left(AM\right)y=b\begin{array}{c}\end{array}\left(\text{right preconditioned},\text{where}x=My\right).\end{array}$

A good preconditioner shall render AM or MA much better conditioned than A (i.e., M should approximate, in a sense,A ^-1) and, in order to be of practical interest, its construction and application must not be overly expensive.

We now present our preconditioner, which combines ideas from domain decomposition methods, DDM, and level-based Incomplete LU factorization, ILU(k). In DDM, one tries to build an approximation to the action of A ^-1 based on the (approximation of) inverses of smaller, local versions of A (subdomain-interior submatrices A _JJ , in our case). In this work, the action of the local inverses is approximated by ILU(k _Int), while other components required by the preconditioner (see equation (3.4)) are approximated with different levels (k _Bord, k _Prod or k _Γ). The motivation for this approach is that ILU(k) has been widely used with success as a preconditioner in reservoir simulation and DDM has been shown to be an efficient and scalable technique for parallel architectures.

It is well established that in linear systems associated with parabolic problems, the high-frequency components of the error are damped quickly, while the low-frequency ones take many iterations to fade, see³¹. In reservoir simulation, the pressure unknown is of parabolic nature. A two-level preconditioner tackles this problem by combining a “fine” preconditioner component M _F , as the one mentioned in the previous paragraph, with a coarse component M _C , the purpose of which is to annihilate the projection of the error onto a coarse space (associated with the low frequencies). If the two components are combined multiplicatively (see²⁹), the resulting two-level preconditioner is

$M=M_F+M_C-M_{F}A{M}_C.$ (3.1)

In Subsection 3.1 we present a ILU(k)-based fine component M _F and in Subsection 3.2 we describe a simple coarse component.

3.1 ILU(k) based Domain Decomposition

The fine component of the domain decomposition preconditioner is based on the following block LU factorization of A:

$A=LU=\left[\begin{array}{cccc}{L}_1& & & \\ & \ddots & & \\ & & L_P& \\ B_1& \cdots & B_P& I\\ & & & \end{array}\right]\left[\begin{array}{cccc}{U}_1& & & C_1\\ & \ddots & & ⋮\\ & & U_P& C_P\\ & & & S\\ & & & \end{array}\right],$ (3.2)

where A _JJ = L _J U _J is the LU factorization of A _JJ , $B_J=A_{\Gamma J}{U}_J^{-1},C_J=L_J^{-1}{A}_{J\Gamma }$ and

$S=A_{\Gamma \Gamma }-{\displaystyle \sum _{J=1}^P{A}_{\Gamma J}{A}_{JJ}^{-1}{A}_{J\Gamma }=A_{\Gamma \Gamma }-{\displaystyle \sum _{J=1}^P{B}_J{C}_J,}}$ (3.3)

is the Schur complement of A with respect to the interior points.

From the decomposition in (3.2), we can show that the inverse of A is

$$\begin{gathered} A^{-1}=\left[\begin{array}{cccc}{U}_1^{-1}& & & \\ & \ddots & & \\ & & U_P^{-1}& \\ & & & I\\ & & & \end{array}\right]\left[\begin{array}{cc}& C_1\\ I& ⋮\\ & C_P\\ & -I\\ & \end{array}\right]\left[\begin{array}{cc}& \\ I& \\ & \\ & S^{-1}\\ & \end{array}\right] \\ \left[\begin{array}{cccc}& & & \\ & I& & \\ & & & \\ B_1& \cdots & B_P& -I\\ & & & \end{array}\right]\left[\begin{array}{cccc}{L}_1^{-1}& & & \\ & \ddots & & \\ & & L_P^{-1}& \\ & & & I\\ & & & \end{array}\right]. \end{gathered}$$ (3.4)

We want to define a preconditioner M _F approximating the action of A ^-1 on a vector. Therefore, we need to define suitable approximations for $L_J^{-1},U_J^{-1}$, B _J and C _J , J=1,…,P, and for S ^-1. These approximations are denoted by ${\tilde{L}}_J^{-1},{\tilde{U}}_J^{-1},{\tilde{B}}_J,{\tilde{C}}_J$, and $S_F^{-1}$, respectively. The fine preconditioner is then defined as

$$\begin{gathered} M_F=\left[\begin{array}{cccc}{\tilde{U}}_1^{-1}& & & \\ & \ddots & & \\ & & {\tilde{U}}_P^{-1}& \\ & & & I\\ & & & \end{array}\right]\left[\begin{array}{cc}& {\tilde{C}}_1\\ I& ⋮\\ & {\tilde{C}}_P\\ & -I\\ & \end{array}\right]\left[\begin{array}{cc}& \\ I& \\ & \\ & S_F^{-1}\\ & \end{array}\right] \\ \left[\begin{array}{cccc}& & & \\ & I& & \\ & & & \\ {\tilde{B}}_1& \cdots & {\tilde{B}}_P& -I\\ & & & \end{array}\right]\left[\begin{array}{cccc}{\tilde{L}}_1^{-1}& & & \\ & \ddots & & \\ & & {\tilde{L}}_P^{-1}& \\ & & & I\\ & & & \end{array}\right]. \end{gathered}$$ (3.5)

In the remaining of this subsection, we describe precisely how these approximations are chosen.

First we define ${\tilde{L}}_J$ and ${\tilde{U}}_j$ as the result of the incomplete LU factorization of A _JJ with level of fill k _Int, $\left[{\tilde{L}}_J,{\tilde{U}}_J\right]=ILU\left(A_{JJ},k_{Int}\right)$. In our numerical experiments, we used k _Int = 1, which is a usual choice in reservoir simulation. Even though ${\tilde{L}}_J$ and ${\tilde{U}}_J$ are sparse, ${\tilde{L}}_J^{-1}{A}_{J\Gamma }$ and ${\tilde{U}}_J^{-T}{A}_{\Gamma J}^T$ (which would approximate C _J and $B_J^T$) are not. To ensure sparsity, ${\tilde{C}}_J\approx {\tilde{L}}_J^{-1}{A}_{J\Gamma }$ is defined as the result of the incomplete triangular solve with multiple right-hand sides of the system ${\tilde{L}}_J{\tilde{C}}_J=A_{J\Gamma }$, by extending the definition of level of fill as follows. Let v _l and w _l be the sparse vectors corresponding to the l-th columns of A _JΓ and $\tilde{C}$ respectively, so that ${\tilde{L}}_J{w}_l=v_l$. Based on the solution of triangular systems by forward substitution, the components of v _l and w _l are related by

$w_{l_k}=v_{l_k}-{\displaystyle \sum _{i=1}^{k-1}{\tilde{L}}_{J_{ki}}{w}_{l_i}.}$ (3.6)

The level of fill-in of component k of w _l is defined recursively as

$Lev\left(w_{l_k}\right)=\min \left\{Lev\left(v_{l_k}\right),\underset{1\le i\le \left(k-1\right)}{\min }\left\{Lev\left({\tilde{L}}_{J_{ki}}\right)+Lev\left(w{l}_i\right)+1\right\}\right\},$ (3.7)

where $Lev\left(v_{l_k}\right)=0$ when $v_{l_k}\ne 0$ and $Lev\left(v_{l_k}\right)=\infty$ otherwise, and $Lev\left({\tilde{L}}_{J_{ki}}\right)$ is the level of fill of entry ki in the ILU(k _Int) decomposition of A _JJ when ${\tilde{L}}_{J_{ki}}\ne 0$ and $Lev\left({\tilde{L}}_{J_{ki}}\right)=\infty$ otherwise. The approximation ${\tilde{C}}_J$ to $C_J=L_J^{-1}{A}_{J\Gamma }$ is then obtained by an incomplete forward substitution (IFS) with level of fill k _Bord, in which we do not compute any terms with level of fill greater than k _Bord during the forward substitution process. We denote ${\tilde{C}}_J=IFS\left({\tilde{L}}_J,A_{J_{\Gamma }},k_{Bord}\right)$. Notice that when k _Bord = k _Int, ${\tilde{C}}_J$ is what would result from a standard partial incomplete factorization of A. The approximation ${\tilde{B}}_J$ to $B_J=A_{\Gamma J}{\tilde{U}}_J^{-1}$ is defined analogously.

Similarly, in order to define an approximation $\tilde{S}$ to S, we start by defining $F_J={\tilde{B}}_J{\tilde{C}}_J$ and a level of fill for the entries of F _J ,

$Lev\left(F_{J_{kl}}\right)=\min \left\{Lev\left(A_{\Gamma {\Gamma }_{kl}}\right),\underset{1\le i\le m}{\min }\left\{Lev\left({\tilde{B}}_{J_{ki}}\right)+Lev\left({\tilde{C}}_{J_{il}}\right)+1\right\}\right\},$ (3.8)

where $m=\#{\Omega }_J^{Int}$ is the number of columns in ${\tilde{B}}_J$ and rows in ${\tilde{C}}_J,Lev\left(A_{\Gamma {\Gamma }_{kl}}\right)=0$ when $A_{\Gamma {\Gamma }_{kl}}\ne 0$ and $Lev\left(A_{\Gamma {\Gamma }_{kl}}\right)=\infty$, otherwise, and $Lev\left({\tilde{C}}_{J_{ki}}\right)$ is the level of fill according to definition (3.7) when ${\tilde{C}}_{J_{ki}}\ne 0$ and $Lev\left(C_{J_{ki}}\right)=\infty$, otherwise ($Lev\left({\tilde{B}}_{J_{il}}\right)$ is defined in the same way). Next, ${\tilde{F}}_J$ is defined as the matrix obtained keeping only the entries in F _J with level less than or equal to k _Prod according to (3.8). We refer to this incomplete product(IP) as ${\tilde{F}}_J=IP\left({\tilde{B}}_J,{\tilde{C}}_J,k_{\mathrm{Pr}od}\right)$.

$\tilde{S}$ is then defined as

$\tilde{S}=A_{\Gamma \Gamma }-{\displaystyle \sum _{J=1}^P{\tilde{F}}_J.}$ (3.9)

While $\tilde{S}$ approximates S, we need to define an approximation for S ^-1. Since $\tilde{S}$ is defined on the global interface Γ, it is not practical to perform ILU on it. Instead, we follow the approach employed in⁷ and define for each subdomain a local version of $\tilde{S}$,

${\tilde{S}}_J=R_J\tilde{S}{R}_J^T,$ (3.10)

where $R_J:\Gamma \to {\overline{\Gamma }}_J$ is a restriction operator such that ${\tilde{S}}_J$ is the result of pruning $\tilde{S}$ so that only the rows and columns associated with ${\overline{\Gamma }}_J$ remain. More precisely, if $\left\{i_1,i_2\mathrm{,...,}{i}_{n_{{\overline{\Gamma }}_J}}\right\}$ is a list of the nodes in Γ that belong to ${\overline{\Gamma }}_J$, then the k-th row of R _J is $e_{i_k}^T$, the i _k -th row of the n _Γ × n _Γ identity matrix,

$R_J=\left[\begin{array}{c}{e}_{i_1}^T\\ ⋮\\ e_{i_{n_{\Gamma }}}^T\end{array}\right].$

Finally, our approximation $S_F^{-1}$ to S ^-1 is defined as

$S_F^{-1}={\displaystyle \sum _{J=1}^P{T}_J{\left(L_{{\tilde{S}}_J}{U}_{{\tilde{S}}_J}\right)}^{-1}{R}_J\approx }{\displaystyle \sum _{J=1}^P{T}_J{\tilde{S}}_J^{-1}{R}_J,}$ (3.11)

where $L_{{\tilde{S}}_J}$ and $U_{{\tilde{S}}_J}$ are given by ILU(k _Γ) of ${\tilde{S}}_J$. Here $T_J:{\overline{\Gamma }}_J\to \Gamma$ is an extension operator that takes values from a vector that lies in ${\overline{\Gamma }}_J$, scales them by $w_1^J\mathrm{,...,}{w}_{n_{{\overline{\Gamma }}_J}}^J$ (which we call weights), and places them in the corresponding position of a vector that lies in Γ. Therefore, using the same notation as before, the k-th column of T _J is $w_k^J{e}_{i_k}$,

$T_J=\left[w_1^J{e}_{i_1}\begin{array}{c}\end{array}\cdots \begin{array}{c}\end{array}{w}_{n_{\Gamma }}^J{e}_{i_{n_{\Gamma }}}\right].$

Different choices for the weights gives rise to different options for T _J , our choice is one that avoids communication among processes and is defined as follows.

$w_k^J=\{\begin{array}{c}\mathrm{1,}\text{if the}{i}_k\text{-th node of}\Gamma \text{belongs to}{\Gamma }_J\begin{array}{cc}& \end{array}\\ \mathrm{0,}\text{if the}{i}_k\text{-th node of}\Gamma \text{belongs to}{\overline{\Gamma }}_J-{\Gamma }_J.\end{array}$

We note that T _J ,J=1,...,P, form a partition of unity. Our approach is based on the Restricted Additive Schwarz method, proposed on⁵.

3.2 Coarse Space Correction

We define a coarse space spanned by the columns of an (n×P) matrix that we call $R_0^T$. The J-th column of $R_0^T$ is associated with the extended subdomain ${\overline{\Omega }}_J$ and its i-th entry is

${\left(R_0^T\right)}_{iJ}=\{\begin{array}{c}\mathrm{0,}\text{if node}i\text{is not in}{\overline{\Omega }}_J\text{and}\\ {\mu }_{{\Omega }_i}^{-1},\text{if node}i\text{is in}{\overline{\Omega }}_J,\end{array}$

where the i-th entry of µ _Ω, denoted ${\mu }_{{\Omega }_i}$, counts how many extended subdomains the i-th node belongs to. Notice that ${\left(R_0^T\right)}_{iJ}=1\begin{array}{c}\end{array}\forall i\in {\Omega }_J^{Int}$ and that the columns of $R_0^T$ form a partition of unity, in the sense that their sum is a vector with all entries equal to 1.

We define M _C by the formula

$M_C=R_0^T{\left(R_{0}A{R}_0^T\right)}^{-1}{R}_0.$ (3.12)

Notice that this definition ensures that M _C A is a projection onto range $\left(R_0^T\right)$ and for A symmetric positive definite this projection is A-orthogonal. Since $R_{0}A{R}_0^T$ is small (P×P), we use exact LU (rather than ILU) when applying its inverse.

Summing up, the complete preconditioner has two components: one related to the complete grid, called M _F in equation (3.5), and another related to the coarse space (3.12). It is possible to subtract a third term that improves the performance of the whole preconditioner, although increasing its computational cost. The combined preconditioners are written down as

$M=M_F+M_C\text{or}$ (3.13)

$M=M_F+M_C-M_{F}A{M}_C.$ (3.14)

This formulation implies that the preconditioners will be applied additively (3.13) or multiplicatively (3.1) and can be interpreted as having two levels, see⁷. In the following sections we address this two-level algebraic ILU(k)-based domain decomposition preconditioner by iSchur (for “incomplete Schur”).

4 PARALLEL IMPLEMENTATION

In this section we discuss the parallel implementation of the preconditioner, considering data locality, communication among processes, and strategies to avoid communication.

We use distributed memory model, so the matrix A is distributed among the processes and, as we use PETSc, the rows corresponding to the elements of each subdomain, including interior and local interfaces, reside in the same process, see Figure 3 for a four domain representation.

Figure 3:
Subdomain and matrix distribution among processes.

4.1 Preconditioner Setup

Algorithm 1 describes the construction of the fine and coarse components of the preconditioner.

Algorithm 1:
Preconditioner Setup

Step 1 of Algorithm 1 is the incomplete LU factorization. Because each subdomain is completely loaded in only one process and there is no connections between subdomain interiors, this step can be done in parallel in each process without communication. We use a native and optimized PETSc routine that computes the incomplete LU factorization.

In Step 2, ${\tilde{C}}_J=IFS\left({\tilde{L}}_J,A_{J\Gamma }\mathrm{,0}\right)$ computes an approximation ${\tilde{C}}_J={\tilde{L}}^{-1}{A}_{J\Gamma }$ through a level zero incomplete triangular solver applied to the columns of A _JΓ . For each column a _i of A _JΓ , we take ℐ _i , the set of indices of the nonzero elements of a _i . Then we solve $c_i\left({\mathcal{I}}_i\right)=\tilde{L}\left({\mathcal{I}}_i,{\mathcal{I}}_i\right)\backslash a_i\left({\mathcal{I}}_i\right)$, which is a small dense triangular system, where c _i is the i-th column of ${\tilde{C}}_J$. See Figure 4. The light gray rows and columns represent the ℐ _i indices applied to the rows and columns of the matrix ${\tilde{L}}_J$ and the dark gray elements are the ones in $\tilde{L}\left({\mathcal{I}}_i,{\mathcal{I}}_i\right)$. We use a PETSc routine to get the submatrices $\tilde{L}\left({\mathcal{I}}_i,{\mathcal{I}}_i\right)$ as a small dense matrix and use our own dense triangular solver to compute c _i (ℐ _i ), which is a small dense vector, and then use another PETSc routine to store it in the correct positions of c _i , which is a column of the sparse matrix ${\tilde{C}}_J$. We observe that the nonzero columns of A _JΓ that take part in these operations, in each subdomain, are related to the extended interface, but are stored locally.

Figure 4:
Incomplete Forward Substitution.

Note that ${\tilde{L}}_J$ and A _JΓ (see Equation (2.5)) are in process J, as can be seen in Figure 3, which means that Step 2 is done locally without any communication.

Step 3 computes ${\tilde{B}}_J$ also through the incomplete forward substitution using the same idea described for Step 2 , but A _ΓJ is not stored entirely in process J. Actually, A _ΓJ has nonzero elements in the rows associated with the extended interface ${\overline{\Gamma }}_J$, therefore the elements of A _ΓJ are distributed among process J and its neighbors. So each process has to communicate only with its own neighbors to compute ${\tilde{B}}_J$. The matrix ${\tilde{B}}_J$ is stored locally in process J.

The communication layout of Step 3 is illustrated in an example for four subdomains in Figure 5, where the arrows show the data transfer flow. In this example, using, for instance, five-point centered finite differences, each subdomain has at most two neighbors, but in different 2-D or 3-D grids and with different discretization schemes, each subdomain can have more neighbors.

Figure 5:
Domain distribution among processes.

The Step 4 computes the Schur complement through an incomplete product as described by equation (3.9). The submatrix A _ΓΓ is formed by the subdomains interfaces and therefore is distributed among the processes. Each process computes locally its contribution $IP\left({\tilde{B}}_J,{\tilde{C}}_J,k_{Prod}\right)$ and subtracts globally from A _ΓΓ to build up $\tilde{S}$, which is also distributed among the processes in the same way A _ΓΓ is. Again, each subdomain communicates with their neighbors following the pattern illustrated in Figure 5, but with arrows pointers reversed, and there is a synchronization barrier to form the global matrix $\tilde{S}$.

In Step 5 and 6, each process takes a copy of the local part of $\tilde{S}$ relative to the extended interface ${\overline{\Gamma }}_J$ as in (3.10), communicating only with neighbors, and makes an incomplete LU factorization. The use of local copies avoids communication during the parallel application of $\tilde{S}$ throughout the linear solver iterations when using the RAS⁵ version of the preconditioner.

In Step 7 the matrix $R_{0}A{R}_0^T$ is initially formed as a parallel matrix, distributed among the processes. Local copies of this matrix is made in all processes and then a LU factorization is done redundantly. The order of this matrix is P, which is in general quite small.

4.2 Preconditioner Application

Algorithm 2 describes the application of the preconditioner M to the vector r, i.e., the computation of z=Mr. We define two additional restriction operators, $R_{{\Omega }_J}:\Omega \to {\Omega }_J^{Int}$ and R _Γ:Ω→ Γ. These restriction operators are needed to operate on the correct entries of the vector, it can be seen as a subset of indices in the implementation.

Algorithm 2:
Preconditioner Application

Step 1 of Algorithm 2 solves a triangular linear system in each process involving the ${\tilde{L}}_J$ factor and the local interior portion of the residual, obtaining z _J . This process is done locally and no communication is necessary because the residual r is distributed among the processes following the same parallel layout of the matrix.

Step 2 applies ${\tilde{B}}_J$ multiplying z _J and subtracting from r _Γ obtaining z _Γ and the communication between the neighbor subdomains is necessary only in r _Γ and z _Γ because ${\tilde{B}}_J$ is stored locally in process J.

The Schur complement application is done in Step 3. Because each process has a local copy of its parts of $\tilde{S}$ factored into $L_{\tilde{S}}$ e $U_{\tilde{S}}$, only the communication in z _Γ between the neighbor subdomains is necessary. As described above, T _J is a diagonal operator, so we can compute it as a vector, which is applied through an element wise product. The result is then summed globally in z _Γ. At this point there is a synchronization barrier among all processes.

In Step 4 we apply the ${\tilde{U}}_J$ factor to $z_J-{\tilde{C}}_J{z}_{\Gamma }$. The matrix-vector product ${\tilde{C}}_J{z}_{\Gamma }$ requires communication among neighbors in z _Γ and then the triangular solver is done locally with no communication.

The vectors z _J and the vector z _Γ are gathered in z _F in Step 5, where z _F is the residual preconditioned by the fine component of the preconditioner.

Step 6 applies the coarse component, completing the computation of the preconditioned residual. Triangular solvers involved are local because there are local copies of L _C and U _C in each process, so no communication is necessary.

5 COMPARISON TESTS

We use two metrics to evaluate the preconditioners: the number of iterations and the solution time. In this section we describe the platform where the tests were run, present the results of the proposed preconditioner, and compare with a native PETSc preconditioner.

We consider five matrices from reservoir simulation problems. They are linear systems arising from applying Newton’s method to solve the nonlinear equations discretizing the system of PDEs that model multiphase flow in porous media. Matrices F8, U26, U31 and Psalt were generated from a Petrobras reservoir simulation tool. The Psalt problem is a synthetic presalt reservoir model. The fifth matrix is the standard SPE10 model¹⁰ with over a million active cells. All the problems are black oil models using fully implicit method and there are four degrees of freedom per grid cell.

Table 1 shows the problems’ size considering the total number of cells and the size of their respective matrices.

Table 1:
Test case problems.

Matrix	Problem Size	Matrix Size
F8	182558	730232
U26	408865	1635460
U31	617459	2469836
Psalt	765620	3062480
SPE10	1094421	4377684

For both preconditioners the problem was solved using native PETSc right preconditioned GMRES method²⁷, with relative residual tolerance of 1e-4, maximum number of iterations of 1000, and restarting at every 30^th iteration. The mesh was partitioned using PT-Scotch⁹. ISchur has fill-in level 1 in the interior of the subdomains and zero in the interfaces, namely k _Int = 1, k _Bord = k _Prod = k _Γ = 0, as this specific configuration presented a better behavior in preliminary Matlab tests.

We compare the proposed two-level preconditioner with PETSc’s native block Jacobi preconditioner (see³) combined (as in equation (3.1)) with the same coarse component presented in this paper. We refer to this combination simply as block Jacobi. We note that block-Jacobi/ILU is the default preconditioner in PETSc and is a common choice of preconditioner in reservoir simulation, which makes it an appropriate benchmark.

The first set of tests was run on an Intel Xeon(R) 64-bit CPU E5-2650 v2 @ 2.60GHz with 16 cores (hyper-threading disabled), 20MB LLC, 64GB RAM (ECC enabled). The operating system is a Debian GNU/Linux distribution running kernel version 3.2.65. The compiler is Intel Parallel Studio XE 2015.3.187 Cluster Edition.

The iSchur preconditioner, both in setup and application phases, requires more computation and communication than block Jacobi, so it has to present a pronounced reduction in number of iterations compared with block Jacobi in order to be competitive. Table 2 shows the number of iterations and the total time, in seconds, for both preconditioners, for all the tested problems, ranging from one to sixteen processes. For each test, the best results are highlighted in gray. We can observe that the iSchur preconditioner does reduce the number of iterations compared to block Jacobi for all problems in any number of processes. For one process (sequential execution), both preconditioners are equal, as there are no interface points, so we show the results just for the sake of measuring speed up and scalability. In average, iSchur does 80% as many iterations as block Jacobi. In this set of experiments, iSchur was more robust than block Jacobi, as the latter failed to converge in two experiments.

Table 2:
Number of iterations and the total time for both preconditioners from 1 to 16 processes. (NC) stands for “not converged”.

Mat	Prec	Number of Iterations Processes					Total Time (s) Processes
		1	2	4	8	16	1	2	4	8	16
F8	iS	104	0.9113	0.9106	0.992	0.9102	10.8	7.9	4.1	2.1	1.8
	bJ		119	109	125	135		0.96.1	0.93.0	0.91.9	0.91.7
U26	iS	92	0.9131	0.9199	0.9218	0.9303	22.0	0.919.8	0.914.5	0.99.1	9.9
	bJ		252	NC	404	349		28.2	NC	13.3	0.99.7
U31	iS	69	0.9110	0.9449	0.9517	0.9864	28.0	28.8	0.948.7	0.932.1	0.943.1
	bJ		124	593	639	NC		0.924.1	56.7	35.4	NC
Psalt	iS	34	0.968	0.957	0.972	0.954	17.8	24.3	11.9	8.2	5.6
	bJ		87	76	117	76		0.919.7	0.99.2	0.98.0	0.95.0
SPE10	iS	116	0.9135	0.9148	0.9148	0.9176	71.5	53.9	30.9	18.0	17.5
	bJ		142	149	164	198		0.943.6	0.924.1	0.915.4	0.916.1

We also observe that, for some of the matrices, the number of iterations increases with the number of processes. This is to be expected, since block Jacobi neglects larger portions of the matrix as the size of the subdomains decreases. A correlated fact occurs with iSchur. In this case as we are using different levels of fill in the interface related parts and in the interior parts, more information is neglected with more subdomains, as the number of interface points increases.

Table 2 also shows the total time of the solution using both preconditioners. We can see that despite the reduction in the number of iterations the solver time of iSchur is smaller than block Jacobi only in a few cases. One of the main reasons is that the iSchur setup time is greater than block Jacobi’s, especially when running on few processes. Although iSchur’s setup is intrinsically more expensive than block Jacobi’s, we believe that this part of our code can still be further optimized to reduce this gap.

In reservoir simulation, the sparsity pattern of the Jacobian matrix only changes when the operating conditions of the wells change, and therefore it remains the same over the course of many nonlinear iterations. This is an opportunity to spare some computational effort in the setup phase, since in this case the symbolic factorization does not need to be redone. We can see in Table 3 the time comparison when the symbolic factorization is not considered; in this scenario iSchur is faster in most cases. Furthermore, block Jacobi is a highly optimized, mature code (the default PETSc preconditioner), while iSchur could possibly be further optimized.

Table 3:
Time of iterative solver without the symbolic factorization.

Iterative Solver Time (s)
Mat	Prec	Processes
		1	2	4	8	16
F8	iS	10.5	6.3	3.3	0.91.7	0.91.6
	bJ		0.96.1	0.93.0	1.9	1.7
U26	iS	21.7	0.916.2	0.912.7	0.98.2	0.99.4
	bJ		28.1	NC	13.3	9.7
U31	iS	27.5	0.922.8	0.945.6	0.930.5	0.942.2
	bJ		23.9	56.7	35.3	NC
Psalt	iS	17.3	0.917.5	0.98.4	0.96.4	0.94.6
	bJ		19.6	9.2	8.0	5.0
SPE10	iS	70.8	44.2	26.0	15.5	0.916.1
	bJ		0.943.4	0.924.0	0.915.3	16.1

In order to examine the behavior of the preconditioners as the number of subdomains increases, we ran tests up to 512 processes on a cluster consisting of nodes with 16GB RAM and two Intel Xeon E5440 Quad-core processors each, totaling 8 CPU cores per node. The results are depicted in Figure 6. We see that iSchur consistently outperforms block Jacobi in terms of number of iterations. Nonetheless, the behaviour of both preconditioners, with respect to scalability, cannot be predicted from these experiments. Tests with more MPI processes, in a larger cluster, are still required.

Figure 6:
Comparison between block Jacobi and iSchur up to 512 processes.

6 CONCLUSIONS AND REMARKS

The presented preconditioner successfully reduced the number of iterations for the target reservoir problems compared to block Jacobi for the full range of partitions tested.

The scalability of both preconditioners is an on-going research as tests in a larger cluster are still necessary, but these initial results indicate a little advantage to iSchur.

For the current implementation the total solution time, setup and iterative solver, of the iSchur is smaller in only a few cases, where one important factor is the inefficiency in the setup step. The block Jacobi in PETSc is highly optimized while our implementation has room for improvement.

Removing the symbolic factorization time, considering the case when it is done once and used multiple times, which is usual in reservoir simulation, the solution using the iSchur preconditioner is faster for most tested cases. It shows that the presented preconditioner can be a better alternative than block Jacobi for the finer component when using two-level preconditioners to solve large-scale reservoir simulation problems. Furthermore, we deem that the application part of the code can still be optimized.

REFERENCES

• ¹
  T. M. Al-Shaalan, H. Klie, A.H. Dogru, & M. F. Wheeler. Studies of robust two stage preconditioners for the solution of fullyimplicit multiphase flow problems. In SPE Reservoir Simulation Symposium, 2-4 February, The Woodlands, Texas,number SPE-118722 in MS, (2009).
• ²
  D. A. Augusto, L. M. Carvalho, P. Goldfeld, I. C. L. Nievinski, J.R.P. Rodrigues, & M. Souza. An algebraic ILU(k) based two-level domain decomposition preconditioner. In Proceeding Series of the Brazilian Society of Computational and Applied Mathematics, 3 (2015), 010093-1-010093-7.
• ³
  S. Balay, S. Abhyankar, M. F. Adams, J. Brown, P. Brune, K. Buschelman, L. Dalcin, V. Eijkhout, W. D. Gropp, D. Kaushik, M. G. Knepley, L. C. McInnes, K. Rupp, B.F. Smith, S. Zampini, H. Zhang, & H. Zhang. PETSc users manual. Technical Report ANL-95/11 - Revision 3.7, Argonne National Laboratory, (2016).
• ⁴
  M. Bollhöfer & V. Mehrmann. Algebraic multilevel methods and sparse approximate inverses. SIAM Journal on Matrix Analysis and Applications, 24(1) (2002), 191-218.
• ⁵
  X. Cai & M. Sarkis. A Restricted Additive Schwarz Preconditioner for General Sparse Linear Systems SIAM Journal on Scientific Computing, 21(2) (1999), 792-797.
• ⁶
  H. Cao, H. Tchelepi, J. Wallis & H. Yardumian. Parallel scalable unstructured CPR-type linear solver for reservoir simulation. In SPE Annual Technical Conference and Exhibition, 9-12 October 2005, Dallas, Texas. Society of Petroleum Engineers, (2005).
• ⁷
  L. M. Carvalho, L. Giraud & P. L. Tallec. Algebraic two-level preconditioners for the Schur complement method, SIAM Journal on Scientific Computing, 22(6) (2001), 1987-2005.
• ⁸
  T. F. Chan, J. Xu & L. Zikatanov. An agglomeration multigrid method for unstructured grids. Contemporary Mathematics, 218 (1998), 67-81.
• ⁹
  C. Chevalier & F. Pellegrini. PT-Scotch: A tool for efficient parallel graph ordering Parallel computing, 34(6) (2008), 318-331.
• ¹⁰
  M. A. Christie & M. J. Blunt. Tenth SPE comparative solution project: A comparison of upscaling techniques in SPE Reservoir Simulation Symposium, (2001).
• ¹¹
  D. A. Collins, J.E. Grabenstetter & P. H. Sammon. A Shared-Memory Parallel Black-Oil Simulator with a Parallel ILU Linear Solver. SPE-79713-MS. In SPE Reservoir Simulation Symposium, Houston, 2003. Society of Petroleum Engineers.
• ¹²
  L. Formaggia & M. Sala. Parallel Computational Fluids Dynamics. Practice and Theory., chapter Algebraic coarse grid operators for domain decomposition based preconditioners. Elsevier, (2002), 119-126.
• ¹³
  S. Gries, K. Stüben, G. L. Brown, D. Chen & D. A. Collinns. Preconditioning for efficiently applying algebraic multigrid in fully implicit reservoir simulations. SPE Journal, 19(04): 726-736, August 2014. SPE-163608-PA.
• ¹⁴
  H. Jackson, M. Taroni & D. Ponting. A two-level variant of additive Schwarz preconditioning for use in reservoir simulation. arXiv preprint arXiv: 1401.7227, (2014).
• ¹⁵
  P. Jenny, S. Lee & H. Tchelepi. Multi-scale finite-volume method for elliptic problems in subsurface flow simulation. Journal of Computational Physics, 187(1) (2003), 47- 67.
• ¹⁶
  P. Jenny, S. Lee & H. Tchelepi. Adaptive fully implicit multi-scale finite-volume method for multi-phase flow and transport in heterogeneous porous media. Journal of Computational Physics, 217(2) (2006), 627-641.
• ¹⁷
  P. Jenny, S. H. Lee & H. A. Tchelepi. Adaptive multiscale finite-volume method for multiphase flow and transport in porous media. Multiscale Modeling & Simulation, 3(1) (2005), 50-64.
• ¹⁸
  A. Manea, J. Sewall & H. Tchelepi. Parallel multiscale linear solver for reservoir simulation. Lecture Notes - slides - 4th SESAAI Annual Meeting Nov 5, 2013, (2013).
• ¹⁹
  A. Muresan & Y. Notay. Analysis of aggregation-based multigrid. SIAM Journal on Scientific Computing, 30(2) (2008), 1082-1103.
• ²⁰
  R. Nabben & C. Vuik. A comparison of deflation and coarse grid correction applied to porous media flow. SIAM J. Numer. Anal., 42 (2004), 1631-1647.
• ²¹
  F. Nataf, H. Xiang & V. Dolean. A two level domain decomposition preconditioner based on local Dirichlet-to-Neumann maps. Comptes Rendus Mathematique, 348(21-22) (2010), 1163-1167.
• ²²
  Y. Notay. Algebraic multigrid and algebraic multilevel methods: a theoretical comparison. Numer. Linear Algebra Appl, 12(5-6) (2005), 419-451.
• ²³
  Y. Notay. Algebraic analysis of two-grid methods: The nonsymmetric case. Numerical Linear Algebra with Applications, 17(1) (2010), 73-96.
• ²⁴
  D. W. Peaceman Fundamentals of Numerical Reservoir Simulation. Elsevier, (1977).
• ²⁵
  A. Quarteroni & A. Valli. Applied and Industrial Mathematics: Venice - 1, 1989, chapter Theory and Application of Steklov-Poincaré Operators For Boundary-Value Problems. Kluwer, (1991), 179-203.
• ²⁶
  Y. Saad. Iterative Methods for Sparse Linear Systems. SIAM, 2nd edition, (2003).
• ²⁷
  Y. Saad & M. H. Schultz. GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems SIAM Journal on scientific and statistical computing, 7(3) (1986), 856-869.
• ²⁸
  Y. Saad & J. Zhang. BILUM: Block versions of multielimination and multilevel ILU preconditioner for general sparse linear systems. SIAM Journal on Scientific Computing, 20(6) (1999), 2103-2121.
• ²⁹
  B. Smith, P. Bjørstad & W. Gropp. Domain Decomposition, Parallel Multilevel Methods for Elliptic Partial Differential Equations. Cambridge University Press, New York, 1st edition, (1996).
• ³⁰
  J. Tang, R. Nabben, C. Vuik & Y. Erlangga. Theoretical and numerical comparison of various projection methods derived from deflation, domain decomposition and multigrid methods. Reports of the Department of Applied Mathematical Analysis 07-04, Delft University of Technology, January 2007.
• ³¹
  U. Trottenberg, C. Oosterlee & A. Schuller. Multigrid. Academic Press, Inc. Orlando, USA, (2001).
• ³²
  C. Vuik & J. Frank. Coarse grid acceleration of a parallel block preconditioner. Future Generation Computer Systems, JUN, 17(8) (2001), 933-940.
• ³³
  J. Wallis, R. Kendall & T. Little. Constrained residual acceleration of conjugate residual methods. In SPE Reservoir Simulation Symposium, number SPE 13563 in 8 th Symposium on Reservoir Simulation, Dallas, 1985. Society of Petroleum Engineers.

• †
  This article is based on work presented at CNMAC2016.

• 1
  Had the condition K>J been dropped from definition (2.3),  would still be a separator set. But it would be unnecessarily large, which would yield a more expensive preconditioner.

• 2
  If A is structurally symmetric, then A ^T and A have the same nonzero structure but are not necessarily equal.

Publication Dates

• Publication in this collection
  Jan-Apr 2018

History

• Received
  28 Nov 2016
• Accepted
  26 Sept 2017