Parallel multilevel restricted Schwarz preconditioners for implicit simulation of subsurface flows with Peng-Robinson equation of state

Let $\Omega\subset\mathbb{R}^{d}$, $d=1,2,3$, be the computational domain on which the PDE system (6) are defined. Then a spatial discretization is performed with a mesh $\Omega_{h}$ of the characteristic size $h>0$. To define the one-level Schwarz preconditioner, we first divide $\Omega_{h}$ into nonoverlapping subdomains $\Omega_{i}$, $i=1,\cdots,N_{p}$, and then expand each $\Omega_{i}$ to $\Omega_{i}^{\delta}$, i.e., $\Omega_{i}\subset\Omega_{i}^{\delta}$, to obtain the overlapping partition. Here, the parameter $\delta>0$ is the size of the overlap defined as the minimum distance between $\partial\Omega_{i}^{\delta}$ and $\partial\Omega_{i}$, in the interior of $\Omega$. For boundary subdomains we simply cut off the part outside $\Omega$. We also denote $H>0$ as the characteristic diameter of $\Omega_{i}$, as shown in the left panel of Figure 1.

Let $N$ and $N_{i}$ denote the number of degrees of freedom associated to $\Omega$ and $\Omega_{j}^{\delta}$, respectively. Let $A\in\mathbb{R}^{N\times N}$ be the Jacobian matrix of the linear system defined on a mesh $\Omega_{h}$

Then we can define the matrices ${R}^{\delta}_{i}$ and ${R}^{0}_{i}$ as the restriction operator from $\Omega_{h}$ to its overlapping and non-overlapping subdomains as follows: Its element $({R}^{\delta}_{i})_{l_{1},l_{2}}$ is either (a) an identity block, if the integer indices $1\leq l_{1}\leq N_{i}$ and $1\leq l_{2}\leq N$ are related to the same mesh point and this mesh point belongs to $\Omega_{j}^{\delta}$, or (b) a zero block, otherwise. The multiplication of ${R}_{i}^{\delta}$ with a $N\times 1$ vector generates a smaller $N_{i}\times 1$ vector by discarding all elements corresponding to mesh points outside $\Omega_{i}^{\delta}$. The matrix ${R}_{i}^{0}\in\mathbb{R}^{N_{i}\times N}$ is defined in a similar way, the only difference to the operator ${R}^{\delta}_{i}$ is that its application to a $N\times 1$ vector also zeros out all of those elements corresponding to mesh points on $\Omega_{i}^{\delta}\backslash\Omega_{i}$. Then the classical one-level additive Schwarz (AS, [29]) preconditioner is defined as

with ${A}_{i}={R}^{\delta}_{i}{A}({R}^{\delta}_{i})^{T}$ and $N_{p}$ is the number of subdomains, which is the same as the number of processors. In addition to that, there are two modified approaches of the one-level additive Schwarz preconditioner that may have some potential advantages for parallel computing. The first version is the left restricted additive Schwarz (left-RAS, [6]) method defined by

and the other modification to the original method is the right restricted additive Schwarz (right-RAS, [4]) preconditioner as follow:

In the above preconditioners, we use a sparse LU factorization or incomplete LU (ILU) factorization method to solve the subdomain linear system corresponding to the matrix ${A}_{i}^{-1}$. In the following, we will denote a Schwarz preconditioner simply by $M^{-1}_{h}$ defined on the mesh $\Omega_{h}$ with the characteristic size $h$, when the distinction is not important.

To improve the scalability and robustness of the one-level additive Schwarz preconditioners, especially when a large number of processors is used, we employ the family of multilevel Schwarz preconditioners. Multilevel Schwarz preconditioners are obtained by combining the single level preconditioner $M^{-1}_{h}$ assigned to each level, as shown in Figure 1. For the description of multilevel Schwarz preconditioners [16, 23, 28], we use the index $j=0,1,...,L-1$ to designate any of the $L\geq 2$ levels. The meshes from coarse to fine are denoted by $\Omega_{h_{j}}$, and the corresponding matrices and vectors are denoted $A_{h_{j}}$ and $X_{h_{j}}$. Let us denote $I_{h_{j}}$ as the identity operator defined on the level $j$, and the restriction operator from the level $j$ to the level $j-1$ be defined by

where $N_{j}$ and $N_{j-1}$ denote the number of degrees of freedom associated to $\Omega_{h_{j}}$ and $\Omega_{h_{j-1}}$. Moreover,

is the interpolation operator from the level $j-1$ to the level $j$.

For the convenience of introduction, we present the construction of the proposed multilevel Schwarz preconditioners from the view of the multigrid (MG) V-cycle algorithm [16, 23, 28, 29]. More precisely speaking, in this sense, at each level $j>0$, the Schwarz preconditioned Richardson method works as the pre-smoother and post-smoother, i.e., $M^{-1}_{h_{j}}$ preconditioning the $\mu_{i}\geq 0$ presmoother iterations and $M^{-1}_{h_{j}}$ preconditioning the $\nu_{i}\geq 0$ postsmoother iterations. In the general multigrid V-cycle framework, let $X^{k}$ be the current solution for the linear system (16), the new solution is computed by the iteration $X^{k+1}=MG(b,L,X^{k})$, as described in Algorithm 1.

The focus of the subsection is on the robustness and efficiency of the proposed algorithm for both standard benchmarks and realistic problems in highly heterogeneous media. Unless otherwise specified, the values of physical parameters used in the test cases are set as follows: $\phi=0.2$, $\mu=11.067\times 10^{-3}\;\mbox{cp}$, $R=8.3147295\;\mbox{J}/(\mbox{mol}\cdot\mbox{K})$, $W=16\times 10^{-3}\;\mbox{kg/mol}$, $T=298\;\mbox{K}$, $p_{c}=4.604\times 10^{6}\;\mbox{pa}$, $T_{c}=190.58\;\mbox{K}$, and $T_{b}=111.67\;\mbox{K}$.

We first present results from a 2D test case (denoted as Case-1) in a horizontal layer with the permeability tensor $\mathbf{K}$ in (1) as follows:

With the media being horizontal, we neglect the effect of the gravity. In the configuration, the distribution of permeabilities includes two domains with different isotropy, as shown in Figure 3. The computational domain is 100 meters long and 100 meters wide. In the simulation, we assume that the left and right boundaries of the domain are impermeable. Then we flood the system by gas from the top to the bottom, i.e., we set $p_{in}=10$ atm at the top boundary and set $p_{out}=1$ atm at the right boundary. And there is no injection/extraction inside the domain. Compared with the previous example, in this test case there is an $H$-shape zone and the value of permeabilities has a huge jump inside and outside the zone, which brings about even greater challenges to the fully implicit solver. In the test, the simulation is performed on a $512\times 512$ mesh, the time step size is fixed to $\Delta t=0.1$ day, and the simulation is stopped at $10$ day. Figure 3 also illustrates the contour plots of the pressure. Table 1 shows the performance of the proposed fully implicit solver with respect to different permeability configurations. It is clearly seen that the simulation spends more computing time when the isotropic medium is used, which attributes to the increase in nonlinearity of the problem that affects the number of linear iterations.

In the following, the experiment is conducted to simulate some 3D test cases. The focus is on the flow model with the medium being highly heterogeneous and isotropic, which significantly increases the nonlinearity of the system and imposes a severe challenge on the fully implicit solver. We first consider a 3D domain with dimension $\Omega=(0,100\;\mbox{m})^{3}$, in which the permeabilities of the porous medium are random, denoted as Case-2. The random distribution of permeability with the range [3.1,14426.2] is generated by a geostatistical model using the open source toolbox MRST [19], as shown in Figure 4. In the test, two different flow configurations are used for the injection boundary, i.e., the fixed injection pressure and the fixed injection rate. In the fixed pressure configuration, the flow with the pressure of 10 atm is injected at the part of the left-hand side $x=0\;\mbox{m},y\in[0,10\;\mbox{m}],z\in[0,10\;\mbox{m}]$; while, for the case of the fixed injection rate, the flow is injected from the same sub-boundary with the Darcy’s velocity $\mathbf{v}=5\times 10^{-4}~{}\mbox{m/s}$. The fluid flows from the right-hand side of the domain with a fixed pressure $p_{out}=1$ atm, and no-flow boundary conditions are imposed on the other boundaries of the domain. The initial pressure of $p_{0}=1$ atm is specified in the whole domain and the parameters required for this example are consistent with the previous examples.

In this case, the domain is highly heterogeneous, leading to the increase in nonlinearity of the problem that is challenging for the numerical techniques. Figure 5 shows the plots of the pressure profiles for Case-2 at different times under the two flow configurations. It is demonstrated from the figures that the proposed approach successfully resolves the different stages of the simulation, and the flow is close to reach the breakthrough at the time $t=5$ days. In Table 2, we present the history of the values of nonlinear and linear iterations and the total execution time for the proposed solver at different times. The simulation is carried out on a $128\times 128\times 128$ mesh and the time step size is $\Delta t=0.1$ day. It can be observed from Table 2 that the number of nonlinear and linear iterations has barely changed, while the total execution time increases constantly with the advance of time as expected, which displays the robustness and effectiveness of the implicit solver to handle the highly variety of physical parameters in the model problem.

In the all of the above test cases, the porosity $\phi$ in the model problem (1) is assumed to a constant. In the next 3D test case (Case-3) of this subsection, we import the porosity and the permeability from the Tenth SPE Comparative Project (SPE10) [10] as an example of a realistic realization with geological and petrophysical properties, in which the porosity and the permeability are capable of variation with the change of the position. It is a classic and challenging benchmark problem for reservoir due to highly heterogeneous permeabilities and porosities. As shown in Figure 6, the permeability is characterized by variations of more than six orders of magnitude and is ranged from $6.65\times 10^{-4}$ to $2\times 10^{4}$, and the porosity scale ranges from 0 to 0.5. The 3D domain dimensions are 1200 ft long $\times$ 2200 ft wide $\times$ 170 ft thick. In the test, the boundary condition for pressure on the left-hand side of the domain is uniformly imposed to $p_{in}=10$ atm and on the opposite side is imposed to $p_{out}=1$ atm. No-flow boundary condition is set to other boundary of the domain. Other used parameters are the same as the previous case.

In Figure 7, we display the distributions for the pressure when the simulation is finished at different time points $t=1,3,6,10$ days. In the simulation, the mesh is $60\times 220\times 85$ and the time step size is fixed to $0.1$ day. The results shown in Figure 7 demonstrate that the proposed approach successfully resolves the rapid and abrupt evolution of the simulation at different stages, while keeping the solver in a robust and efficient way. Finally, we analyze the behavior of the proposed fully-implicit method when the time step size $\Delta t$ is changed. In the test, we again run the SPE-10 model on a fixed $60\times 220\times 85$ mesh. The simulation is stopped at $t=3$ year. The results on the average numbers of nonlinear and linear iterations as well as the total compute time are summarized in Table 3. The results in the table clearly indicate that the combination of nonlinear and linear iterations works well for even very large value of time steps. The implicit approach converges for all time steps and is unconditionally stable. In addition to that, we also notice that, as the time step size $\Delta t$ decreases, the average number of nonlinear and linear iterations become smaller, whereas the total computing time increases. This behavior is somehow expected for the fully implicit approach in a variety of applications [2, 16, 32, 35, 36].

In this subsection, we focus on the parallel performance of the proposed fully implicit solver with respect to the one-level or multilevel Schwarz preconditioner by Case-1, Case-2, and a new 3D test case. In the new 3D case, denoted as Case-4, the computational domain is $\Omega=(0,100\;\mbox{m})^{3}$ with the medium being heterogeneous and isotropic, and the distribution of permeability in each layer $z\in[25\;\mbox{m},75\;\mbox{m}]$ is the same as the distribution of Case-1. We flood the system by gas from the behind face ($x\in[0,100\;\mbox{m}]$, $y=100\;\mbox{m}$, $z\in[0,100\;\mbox{m}]$) to the front face, i.e., we set $p_{in}=10$ atm at the behind face and set $p_{out}=1$ atm at the front face. And there is no injection/extraction inside the domain.

In the test, we use the following stopping parameters and notations, unless specified otherwise. The relative and absolute tolerances in (15) for Newton iterations are set to $10^{-6}$ and $10^{-10}$, respectively. The linear systems are solved by the one-level or two-level Schwarz preconditioned GMRES method with absolute and relative tolerances of $10^{-8}$ and $10^{-5}$ in (14), except for the coarse solve of the two-level preconditioner, for which we use $2\times 10^{-1}$ as the relative stopping tolerance. The restart value of the GMRES method is fixed as 30. In the tables of the following tests, the symbol “$N_{p}$” denotes the number of processor cores that is the same as the number of subdomains, “N. It” stands for the average number of Newton iterations per time step, “L. It” the average number of the one-level or two-level Schwarz preconditioned GMRES iterations per Newton iteration, and “Time” the total computing time in seconds.