The multiscale perturbation method for second order elliptic equations

Abstract

In the numerical solution of elliptic equations, multiscale methods typically involve two steps: the solution of families of local solutions or multiscale basis functions (an embarrassingly parallel task) associated with subdomains of a domain decomposition of the original domain, followed by the solution of a global problem. In the solution of multiphase flow problems approximated by an operator splitting method one has to solve an elliptic equation every time step of a simulation, that would require that all multiscale basis functions be recomputed.

In this work, we focus on the development of a novel method that replaces a full update of local solutions by reusing multiscale basis functions that are computed at an earlier time of a simulation. The procedure is based on classical perturbation theory. It can take advantage of both an offline stage (where multiscale basis functions are computed at the initial time of a simulation) as well as of a good initial guess for velocity and pressure. The formulation of the method is carefully explained and several numerical studies are presented and discussed. They provide an indication that the proposed procedure can be of value in speeding-up the solution of multiphase flow problems by multiscale methods.

Introduction

Problems of interest to the oil industry, such as the petroleum reservoirs in Brazil’s pre-salt layer, may require spatial discretizations involving billions of cells such that numerical solutions of multiphase flows can accurately reflect the heterogeneous nature of the underlying subsurface formations. Frequently used reservoir simulators cannot handle such large computational tasks. Multiscale methods are a class of numerical procedures that have good potential to be used in producing approximations for such challenging problems. These methods have been developed both in the context of finite volume and finite element approximations [1], [2]. See [3], [4] for references for more recent work. In the context of mixed finite elements we mention [5], [6], [7] that have been generalized by the recently introduced Multiscale Robin Coupled Method [3], [4]. The aforementioned multiscale mixed methods share a similar algorithmic structure. Initially a set of local solutions (or mixed multiscale basis functions, MMBF) of boundary value problems (BVPs, the methods differ primarily in the choice of boundary conditions for these local problems) is solved for non-overlapping subdomains of the domain of interest. Then, a global problem is defined, that will couple the multiscale basis functions to produce the final approximate solution for the global problem at hand.

We are concerned with the development of the Multiscale Perturbation Method (MPM), a novel procedure to speed-up the parallel (multi-core) numerical solution of multiphase flow problems by multiscale methods that share the structure just described. In order to motivate the discussion of the proposed procedure for second order elliptic equations, we consider the numerical approximation of two-phase (water-oil) flows in porous media by an operator splitting method. Operator splitting methods [8], [9], [10], [11] are well known strategies to handle the coupled systems of partial differential equation by decomposing them into smaller problems. The problems that result from this decomposition describe simpler physics and numerically they are solved sequentially. The governing system of equations for two phase flows is given by (see [9], [12])

$$-\lambda \left(s_w\right)\mathbf{\kappa }\nabla p=\mathbf{u},$$

$$\nabla ·\mathbf{u}=q,$$

$$\varphi \frac{\partial s_w}{\partial t}+\nabla ·\left(\frac{{\lambda }_w\left(s_w\right)}{\lambda \left(s_w\right)}\mathbf{u}\right)=0,$$

where the unknowns are the velocity u(x, t) and pressure p(x, t) fields, and the water saturation s_w(x, t), as a function of space $\mathbf{x}\in \Omega \subset {\mathbb{R}}^d,$ $d=2,3$ and time t ≥ 0. For simplicity in the presentation additional assumptions considered in this model are the absence of capillary pressure and gravity effects, along with a fully saturated media; ($s_w+s_o=1$). The procedure introduced here could incorporate both gravity effects and capillary pressure, in that one would still need to solve a second order elliptic equation for the global pressure (see [9], [12]). The assumption of a fully saturated media is typical of oil-water flows in petroleum reservoirs. In the above equations, $q=q(\mathbf{x},t)$ is a volumetric source term for wells, $\varphi =\varphi \left(\mathbf{x}\right)$ is the porosity, κ(x) is the permeability tensor (we assume it is a scalar). The coefficient λ(s_w) is the total phase mobility, given by

$$\lambda \left(s_w\right)=\frac{{\mathbf{\kappa }}_{rw}\left(s_w\right)}{{\mu }_w}+\frac{{\mathbf{\kappa }}_{ro}\left(s_w\right)}{{\mu }_o},$$

where ${\mathbf{\kappa }}_{r\alpha }\left(s_w\right),(\alpha =w,o)$ are the relative permeabilities and μ_α is the viscosity of the α phase. We define the α phase mobility by

$${\lambda }_{\alpha }\left(s_w\right)=\frac{{\mathbf{\kappa }}_{r\alpha }\left(s_w\right)}{{\mu }_{\alpha }},$$

such that $\lambda \left(s_w\right)={\lambda }_o\left(s_w\right)+{\lambda }_w\left(s_w\right)$. We define $K(\mathbf{x},s_w)=\lambda \left(s_w\right)\mathbf{\kappa }\left(\mathbf{x}\right)$. Initial conditions for one of the saturations and boundary conditions have to be provided.

For computational efficiency an operator splitting method is frequently used for the numerical approximation of the above system. This splitting is performed by “freezing” the saturation in time for the solution of the second order elliptic equation, allowing for the segregated solution of both equations. In one time step of the overall method, say from time t_n to time $t_{n+1},$ saturation is kept at time t_n for the solution of Eqs. (1) and (2), to obtain $p^{n+1}$ and ${\mathbf{u}}^{n+1}$ both at time $t_{n+1}$. These newly computed pressure and velocity are then used in Eq. (3) to evolve saturation to the next time step $t_{n+1}$. The main advantage is that now each equation can be solved by well known techniques that can benefit from different time steps. A larger time step is used for the update of velocity and pressure and smaller time steps, that satisfy a CFL condition [13], are used in the solution of a hyperbolic conservation law for the water saturation. The use of distinct time steps in illustrated in Fig. 1.

The solution of two phase flows can naturally take advantage of multiscale methods that have been developed for second order elliptic equations by replacing a global solver by a multiscale one in the update of velocity and pressure. The proposed procedure aims at finding the solution for velocity and pressure by replacing the construction of a set of multiscale basis functions or local solutions by reusing the multiscale basis functions computed for another elliptic equation. Thus, the proposed method is formulated in terms of two elliptic boundary value problems (let us refer to them as the initial time and target problems) that are approximated in terms of a common decomposition of the domain. For the initial time problem a set of multiscale basis functions is computed (say, in the first elliptic solution of a time dependent problem). Within the MPM the solution of the target elliptic equation is obtained in two steps. In the first step an interface problem is solved: perturbation theory is applied to predict the normal component of the fluxes on the skeleton of the domain decomposition. Through perturbation theory we derive a hierarchy of BVPs that share the same coefficients with the time zero problem. Thus, these BVPs can be solved by taking advantage of the pre-computed, time zero MMBFs. The first step is followed by a task that fits well in multi-core devices: for each subdomain a BVP is solved using the boundary conditions computed in the first step, giving the final approximation for the velocity field.

The proposed method has some distinctive properties that can be summarized as follows:

• The method works quite well for high-contrast formations;

• A considerable fraction of the computational work can be performed in an embarrassingly parallel offline stage (all solutions of local problems can be naturally solved independently from each other);

• If a full set of multiscale basis functions (one multiscale basis function for each edge of the boundary of a subdomain) is computed for the time zero problem the target solution has continuous fluxes (thus, a downscaling step can be avoided);

• The method can incorporate a good initial guess (that would be relevant in the solution of time dependent problems);

• The cost of the procedure is considerably smaller than a typical multiscale method;

• The procedure in non-intrusive: it can incorporate various multiscale methods in their original formulation and numerical implementation;

• Numerical studies indicate that the accuracy of the procedure is comparable (or better) to that of the multiscale method used in its formulation.

This work is organized as follows. In Section 2 we recall the Multiscale Mixed Method [7] that will be the underlying multiscale method for the description of the proposed procedure. Then, the new method is derived in detail in Section 3. Initially we discuss a mixed (for the flux and pressure variables) perturbation method for a BVP. Then we derive an auxiliary residual problem that allows one to produce improved solutions by taking advantage of an initial guess (that is frequently available in the solution of time dependent methods). Section 4 is dedicated to a discussion of numerical solutions of a model problem. Our conclusions and the prospects of future work appear in Section 5. The dimensionless form of the Poisson problem appears in an appendix.