Multiscale Mixed Methods with Improved Accuracy: The Role of Oversampling and Smoothing

Multiphase flows in subsurface formations are governed by coupled systems of partial differential equations. In the case of two and three-phase flows of interest to the oil industry frequently these equations are decoupled by operator splitting methods (see [1, 2] for such methods developed, respectively, for two and three-phase flows) and a convection-diffusion equation needs to be solved along with a second-order elliptic equation. The focus of this work is on the numerical solution of these elliptic problems that are challenging because the permeability coefficient may exhibit very large contrast and classical iterative schemes are known to converge slowly.

In this paper, we direct our attention to multiscale methods, which have received significant attention due to their natural parallelizability on multi-core computers. Over the past two decades, various research groups have pioneered specific methodologies integrating domain decomposition with multiscale techniques. These methods generally fall into three distinct categories: Multiscale Finite Volume, Multiscale Finite Element, and Multiscale Mixed Finite Element approaches.

The Multiscale Finite Volume Method was introduced in 1997 to solve elliptic problems encountered in two-phase flows [3, 4]. Following this, the Multiscale Finite Element Method was developed in 2003 [5]. In this study, we focus on Multiscale Mixed Finite Element Methods, which encompass several variants. Here we mention the ones that are directly relevant to our work. The Multiscale Mixed Finite Element Method (MsMFEM) was first formulated in 2003 [6], followed by the development of the Multiscale Mortar Mixed Finite Element Method (MMMFEM) in 2007 [7], and the Multiscale Hybrid-Mixed Finite Element Method (MHM) in 2013 [8]. Our research focuses particularly on the Multiscale Robin Coupled Method (MRCM), introduced in 2018 [9], which can be considered a generalization of the Multiscale Mixed Method (MuMM) introduced in 2014 [10]. This generalization allows for arbitrary and independent interface spaces for the pressure and flux variables.

The exceptional scalability of MRCM for solving industry-relevant subsurface flow problems on multi-core supercomputers has been demonstrated in previous studies [11, 12]. Recent advancements have introduced new interface spaces utilizing the Singular Value Decomposition [13], as well as approaches that take advantage of the underlying physics of problems [14]. While the pressure and normal flux components in MRCM may have discontinuities in the resulting fine grid solution, a post-processing method has been devised to address this issue, ensuring continuous normal flux components [15]. With this enhancement, MRCM becomes applicable for approximating multiphase flows in porous media. Additionally, the approximation of two-phase flows with MRCM has been explored in recent publications, which have introduced accurate and computationally efficient methods based on operator splitting and implicit time-stepping methodologies [16, 17, 18, 19].

Robin-type boundary conditions are imposed by MRCM on all local problems defined on a non-overlapping partition of the domain. Based in a domain decomposition methodology [20], MRCM ensures compatibility across the interfaces of adjacent subdomains through weak continuity of pressure and normal components of the flux. At the heart of the MRCM framework lies a critical parameter, denoted as $\alpha$, which governs the relative importance of the pressure and the normal component of the flux in the specification of Robin boundary conditions. This parameterization enables MRCM to replicate outcomes from both the Multiscale Mortar Mixed Finite Element Method (MMMFEM) and the Multiscale Hybrid-Mixed Finite Element Method (MHM). In the case of MMMFEM, pressure continuity along the interface is preserved at the fine-grid scale, while flux continuity is weakly satisfied on a larger scale, typically corresponding to the size of each subdomain. By assigning a higher weight to pressure continuity via $\alpha$, MRCM converges towards MMMFEM. Conversely, in MHM, flux continuity along the interface is achieved at the fine-grid scale. Adjusting $\alpha$ to prioritize flux continuity leads MRCM to converge towards MHM [9].

However, even with the latest version of MRCM utilizing non-overlapping partitions, a notable resonance error persists along the interface of solutions obtained with multiscale mixed methods. To mitigate this error, the oversampling method has been introduced. Initially employed in the Multiscale Finite Volume Method in 1997 [21] and subsequently substantiated in 1999 [22, 23], the oversampling method has undergone further development in recent years. Additional references concerning recent advancements in Multiscale Finite Volume methods with oversampling are available in [24, 25]. Numerous references exist for Multiscale Finite Element methods incorporating oversampling. References dedicated to elliptic equations can be located in [26, 27, 28, 29, 30, 31, 32, 33] and papers focusing on two-phase flow can be found in [34, 35, 36]. In the area of Multiscale Finite Element methods with oversampling, specialized methodologies such as Extended Multiscale Finite Element [37, 38, 39], Generalized Multiscale Finite Element [40], and Multiscale DG methods [41] have been also explored.

In the realm of multiscale mixed methods, the first application of the oversampling technique in MsMFEM dates back to 2003 [6]. In this work the authors considered rapidly oscillating coefficients to address significant resonance errors along the interface. Subsequent advancements have further refined oversampling techniques to enhance accuracy and performance. In 2012, an innovative approach introduced a special term into the Multiscale Basis Function, using Green’s function to mitigate resonance errors [42]. Expanding on these efforts, a paper in 2017 integrated upscaling techniques into oversampling methodologies, to approximate two-phase flow problems rather than solely focusing on elliptic equations [43]. This approach shows the importance of considering the entire domain in addressing multiscale challenges. Meanwhile, oversampling for the Mixed Generalized Multiscale Finite Element Method was introduced in 2015 [44], with further advancements in 2017 [45]. These works utilized both offline and online stages in the numerical solution to enhance accuracy. The offline stage involves compiling boundary conditions for each multiscale basis function using preset spaces such as piecewise constant or piecewise linear functions. In contrast, the online stage utilizes results from the offline stage to update boundary values.

In our current study, our primary focus lies on enhancing the Multiscale Robin Coupled Method (MRCM) by integrating two distinct strategies aimed at improving the accuracy of numerical solutions while preserving the computational efficiency inherent to MRCM [11, 12]. Our first significant contribution involves the integration of the oversampling technique into MRCM. This novel approach necessitates non-trivial modifications to the original method, as multiscale basis functions now need to be computed in terms of novel informed spaces. Despite these modifications, compatibility conditions are maintained and enforced on the skeleton of an underlying non-overlapping partition of the domain within a nonconforming multiscale approach. Our second contribution involves the introduction of a smoothing step in the framework of multiscale mixed methods, a concept commonly employed in the context of overlapping Schwarz domain decomposition methods [46, 47]. This step serves as a tool to rectify small-scale errors inherent in the multiscale solution, further enhancing the accuracy of our numerical approach. The computational cost of these steps is very small because a factorization computed in the first step can be reused.

Following the description of the novel procedure, we conduct testing in two distinct examples. The first example has an analytical solution, while the second example addresses problems relevant to the oil industry. Through comprehensive studies, we evaluate various aspects including convergence rates, the significance of oversampling, the impact of smoothing steps, the combined effect of oversampling and smoothing, and the influence of the $\alpha$ parameter on numerical solutions. Our findings reveal that in scenarios characterized by modestly sized overlapping regions and utilizing only a few smoothing steps, our new approach yields a remarkable enhancement in flux accuracy, with improvements reaching two orders of magnitude compared to previous methods. Moreover, results produced with the new approach with piecewise-constant interface spaces are comparable to MRCM results obtained with piecewise linear spaces.

The paper is structured into several sections as follows. The first section focuses on the formulation of the new MRCM with oversampling with subsections dedicated to detailing the formulation of the method along with its well-posedness, and the definition of multiscale basis functions. Following this, the next section describes the concept of a smoothing step, outlining its significance and implementation within the context of the MRCM methodology. In the subsequent section, we present the numerical strategy employed to maximize the computational efficiency of the previously introduced methods. The subsequent section is dedicated to numerical studies and is subdivided into subsections that present the results of our numerical experiments. We consider a problem with an analytical solution and another problem involving a permeability field derived from the SPE 10 project. Through these numerical studies, we aim to assess the effectiveness and performance of the proposed enhancements to the MRCM in various scenarios.

We consider single-phase flow in porous media. The governing equations for pressure $p$ and Darcy velocity ${\bf u}$ are given by

where $\Omega\subset\mathbb{R}^{d}$, $d=2$ or $3$ is the domain of the problem, $K$ is a symmetric, uniformly positive definite tensor with components in $L^{\infty}(\Omega)$, $f\,\in\,L^{2}(\Omega)$ is the source term, $g\,\in\,H^{\frac{1}{2}}(\partial\Omega_{p})$ is the pressure condition on the boundary, $z\,\in\,H^{-\frac{1}{2}}(\partial\Omega_{u})$ is the normal velocity condition on the boundary and ${\bf n}$ is the outer normal to $\partial{\Omega}$.

We will define the smoothing steps through local updates of a solution for Eqs. (1 - 4) that has been obtained by the MRCM with oversampling. Let us refer to it as $({\bf u}_{h}^{i,0},p_{h}^{i,0})$, $i=1,\ldots,m$, where the superscript $0$ indicates that this solution is the first one of a sequence that will follow. The steps that produces the solution $({\bf u}_{h}^{i,k},p_{h}^{i,k})$, $i=1,\ldots,m$, $k=1,\ldots,N_{s}$ can be described as follows:

1. 1.

   Set a coloring scheme for oversampling regions within the subdomain partitioning process. The aim is to assign a distinct color to each subdomain to ensure that subdomains sharing the same color do not have any common boundary points. This approach is illustrated in Fig.(2) for a two-dimensional partition (see [47]).

2. 2.

   For each color, loop of all corresponding oversampling regions.

3. 3.

   For all oversampling regions with the same color set the Robin boundary conditions as $\lambda_{h}^{i,k-1}=-\beta_{i}\,{\bf u}_{h}^{i,k-1}\cdot{\bf\check{n}}^{i}|_{\partial\hat{\Omega}_{i}\setminus\partial\Omega}+\pi^{i,k-1}|_{\partial\hat{\Omega}_{i}\setminus\partial\Omega}$ on $\partial\hat{\Omega}_{i}\setminus\partial\Omega$ (see [50]) in terms of existing $({\bf u}_{h}^{i,k-1},p_{h}^{i,k-1})$, $i=1,\ldots,m$ with $\alpha=1$ (other values for $\alpha$ can be used but our numerical studies indicate that this is a good choice).

4. 4.

   For all oversampling regions with the same color simultaneously find $(\hat{\bf u}_{h}^{i},\hat{p}_{h}^{i})$ by solving

   	$\displaystyle(K^{-1}\hat{\bf u}_{h}^{i},{\bf v})_{\hat{\Omega}_{i}}-(\hat{p}_{h}^{i},\nabla\cdot{\bf v})_{\hat{\Omega}_{i}}+(\beta_{i}\,\hat{\bf u}_{h}^{i}\cdot\check{\bf n}^{i},{\bf v}\cdot\check{\bf n}^{i})_{\partial\hat{\Omega}_{i}\setminus\partial\Omega}$
   	$\displaystyle=-(\lambda_{h}^{i,k-1},{\bf v}\cdot\check{\bf n}^{i})_{\partial\hat{\Omega}_{i}\setminus\partial\Omega}-(g,{\bf v}\cdot\check{\bf n}^{i})_{\partial\hat{\Omega}_{i}\cap\partial\Omega_{p}},$		(35)
   	$\displaystyle(q,\nabla\cdot\hat{\bf u}_{h}^{i})_{\hat{\Omega}_{i}}=(f,q)_{\hat{\Omega}_{i}},$		(36)

   hold for all $({\bf v},q)\,\in\,{\bf V}_{h0}^{i}$.

5. 5.

   For all oversampling regions with the same color update $({\bf u}_{h}^{i,k-1},p_{h}^{i,k-1})$ in $\Omega_{i}$ to be the restriction of $(\hat{\bf u}_{h}^{i},\hat{p}_{h}^{i})$ computed above to $\Omega_{i}$.

6. 6.

   After the loop over all colors is completed set $({\bf u}_{h}^{i,k},p_{h}^{i,k})$ to be the existing $({\bf u}_{h}^{i,k-1},p_{h}^{i,k-1})$, $i=1,\ldots,m$.

7. 7.

   Repeat the above steps $N_{s}$ times.

The novel method introduced in this work, comprising MRCM-O followed by smoothing steps, is referred to as the Multiscale Robin Coupled Method with Oversampling and Smoothing (MRCM-OS).

This section provides the numerical strategy employed within our in-house parallel code developed in C/MPI. Algorithm (1) provides a pseudo-code describing the adopted strategy.

Efficiency within our strategy is primarily achieved through the reuse of LDL^T factorizations. Here, we compute the LDL^T-factorization to solve linear problems of the form $Ax=b$ posed on subdomain $\hat{\Omega}_{i}$. It is worth noting that while solving local problems within each subdomain, the matrix $A$ remains constant, with variations occurring solely in vector $b$. As a result, we can compute the LDL^T-factorization of $A$ once and reuse it across all $Ax=b$ calculations. This approach dramatically reduces the computational cost associated with computing the MBFs in Step 2 taking advantage of Step 1 and the subsequent smoothing step in Step 7 taking advantage of Step 6. As previously mentioned, our numerical experiments (not explored in this manuscript) indicate that the choice of $\alpha=1$ for the smoothing leads to improved accuracy with a reduced number of iterations. Additionally, the subsequent numerical investigations in the following section indicate that the smallest errors occur when $\alpha$ is approximately 1. In practice, choosing $\alpha\sim 1$ for MBFs and smoothing can eliminate the need to compute Step 6, allowing the use of the factorization obtained in Step 1 to compute Step 7. This approach enhances accuracy and reduces computational cost. It is noteworthy that the majority of steps are computed concurrently within each MPI-process, underscoring the parallel nature of our approach. The sole exception arises during the computation of the solution to the interface problem in Step 5.

The simulations discussed were executed in an HPC cluster. We discuss the numerical solutions of two model problems: one having a constant permeability coefficient $K\left(\bf x\right)$ with an analytical solution and the other one encompassing a very heterogeneous permeability field from the SPE10 project (http://www.spe.org/web/csp).

For the problem with constant $K\left(\bf x\right)$, we analyze the system of Eq. (1)-(4) subject to trivial Newmann boundary conditions. The permeability field is set to $K\left(\bf x\right)=1$. The source term $q$ is written as $q=8\pi^{2}cos(2\pi x)cos(2\pi y)$. The analytical solution for this problem is given by $p=cos(2\pi x)cos(2\pi y)$. We are going to discuss several partitions of the domain $[0,1]\times[0,1]$. In all these partitions, each subdomain will be discretized by a computational local grid of size $20\times 20$.

For the problem with a heterogeneous permeability field, we consider a rectangular domain and prescribe Dirichlet-type boundary conditions by fixing the pressure to 1 (0) on the left (right) boundary. We set Neumann-type conditions to zero at the top and bottom boundaries. The source term is set as $q(x)=0$. Additionally, the permeability field $K\left(\bf x\right)$ is set to be the 40th (two-dimensional) layer of the SPE10 project, as depicted in Fig.(3). This particular field poses a considerable challenge for numerical solvers due to its variability and the presence of a highly permeable channel. For all reported results in heterogeneous problems, the domain is decomposed into $11\times 3$ subdomains with $H=1/3$, each containing a local $20\times 20$ grid.

We are going to refer to the two problems described above as homogeneous and heterogenous, respectively.

Table 1 summarizes the notation to be considered in the numerical experiments. We remark that the original MRCM is equipped with linear polynomial interface spaces in all studies reported here.

Section 4.1 provides numerical results for both the heterogeneous and homogeneous problems, illustrating some of the advantageous aspects of the proposed MRCM-OS method. First we illustrate our method’s improved accuracy followed by a study showing that MRCM-OS employing only piecewise constant spaces is comparable to the original MRCM with linear polynomial spaces, although with lower computational cost. We emphasize that, in some results of the homogeneous problem, our method utilizing linear polynomial spaces demonstrates the capability to match or even surpass the accuracy of the fine grid solution. Section 4.2 explores the advantages brought by oversampling alone. In Section 4.3, we focus in a detailed analysis of the error reduction achieved by our method employing both oversampling and smoothing techniques.

In conclusion, our study introduces two novel methodologies, oversampling and smoothing, tailored to augment the original MRCM [9], aimed at refining the approximation of subsurface flows. The oversampling approach strategically computes the multiscale basis function within an extended region, effectively mitigating small-scale errors along the interface $\Gamma$ in non-overlapping partitions. Further refinement is achieved through smoothing techniques. Through numerical experiments, our findings demonstrate that our proposed methodology yield a notable improvement in accuracy. Specifically, our method exhibits a two-orders-of-magnitude improvement in flux variables and a one-order-of-magnitude enhancement in pressure, as compared to the original MRCM, across both analytical and heterogeneous problems. Importantly, these improvements are achieved without significant additional computational overhead.

In addition to the aforementioned techniques, our ongoing research has led to the development of novel informed spaces, which shows promise in further minimizing errors, particularly in three-dimensional scenarios. Unlike higher-order polynomial spaces, the informed space offers simplicity in construction while maintaining effectiveness in error reduction. Our forthcoming paper will delve deeper into the exploration and refinement of these spaces, shedding light on its potential contributions to the field of subsurface flow approximation.

The authors acknowledge the National Laboratory for Scientific Computing (LNCC/MCTI, Brazil) for providing HPC resources of the S. Dumont supercomputer, which have contributed to the research results reported within this paper (URL: http://sdumont.lncc.br). Additionally, the computing resources of the Cyber-Infrastructure Research Services at the University of Texas at Dallas Office of Information Technology were utilized. The authors would like to express their gratitude to Dr. Márcio Borges for his assistance in accessing the LNCC cluster.

During the preparation of this work the authors used chatGPT in order to improve language and readability. After using this tool, the authors reviewed and edited the content as needed and take full responsibility for the content of the publication.