Domain decomposition strategies for nonlinear flow problems in porous media

Abstract

Domain decomposition (DD) methods, such as the additive Schwarz method, are almost exclusively applied to linearized equations. In the context of nonlinear problems, these linear systems appear as part of a Newton iteration. However, applying DD methods directly to the original nonlinear problem has some attractive features, most notably that the Newton iterations now solve local problems, and thus are expected to be simpler. Furthermore, strong, local nonlinearities may to a less extent affect the numerical algorithm. For linear problems, DD can be applied both as an iterative solver or as a preconditioner. For nonlinear problems, it has until recently only been understood how to use DD as a solver.

This article offers a systematic study of domain decomposition strategies in the context of nonlinear porous-medium flow problems. The study thus compares four different approaches, which represents DD applied both as a solver and preconditioner, to both the linearized and nonlinear equations. Our model equations are those obtained from a fully implicit discretization of immiscible two-phase flow in heterogeneous porous media. In particular we emphasize the case of nonlinear preconditioning, an algorithm that to our knowledge so far has not been studied nor implemented for flow in porous media. Our results show that the novel algorithm is up to 75% faster than the standard algorithm for the most challenging problems for a moderate number of subdomains.

Introduction

Simulation of flow in subsurface porous media is of crucial importance to a variety of applications, such as oil and gas recovery, CO₂ storage and geothermal energy extraction. In order to give an accurate description of the flow pattern, highly nonlinear partial differential equations are needed. Both the parameter fields (and foremost permeability) and the solution may vary over several orders of magnitudes, and the variations may have different correlation lengths. To resolve all the fine scale variations, a fine discretization is required, which leads to a large number of unknowns. This leads us to seek advanced numerical methods that can handle both high degrees of nonlinearity and multiscale features, for large systems of equations.

The primary variables in reservoir simulation are usually chosen to be the pressure in one phase, and a set of phase saturations or mass variables. The governing equations are an (almost) elliptic pressure equation, and transport equations that are (depending on scale) hyperbolic or parabolic. For realistic physical processes, both the transport and the pressure equations are nonlinear. In order to simulate the spatially discretized nonlinear system without time-step size being limited by stability, an implicit temporal discretization is frequently employed method for some or all equations. With a fully implicit method (FIM), a large nonlinear system must be solved at each time step.

The standard way to handle nonlinearity is by some variation of Newton iterations, in which a system of linear equations must be solved for each iteration. This consumes a considerable amount of the total simulation time, and methods have been developed both to reduce the time it takes to solve linear system, and reduce the number of Newton iterations needed. This involves amongst others two-stage preconditioners [1], [2], [3], reordering techniques [4], [5], and multiscale methods, see e.g. [6], [7], [8].

Another efficient solution strategy for linear and nonlinear equations is domain decomposition (DD) [9], [10], [11], which may be used in conjunction with Newton iterations to enhance the performance of the simulation. This paper is devoted to investigation of how to best apply DD to simulate flow and transport in porous media. Domain decomposition has advantages both for elliptic problems as well as hyperbolic equations [12], and it is therefore well suited for our problems. Importantly, the decoupling of the problem into subdomains makes domain decomposition attractive for parallelization, making efficient large scale computations possible. For nonlinear systems the standard application of domain decomposition techniques is as a linear preconditioner for the solver of the linear system obtained inside the Newton iteration. Conceptually it is also possible to do DD directly to the nonlinear system. The benefit of nonlinear DD is evident: If the nonlinearities mainly are confined to a small region of the computational domain, a spatial decomposition of the nonlinear problem will make the number of Newton iterations on any sub-domain automatically adapt to the degree of nonlinearity. The computational savings of this approach compared to performing Newton iterations on the whole domain in order to reduce the residual in a small region can be substantial. In reservoir simulation, nonlinearities are often mostly confined to advancing fronts; hence nonlinear DD should be of interest to these problems. Common practice has nevertheless been to linearize the system prior to domain decomposition. A major reason for this is that domain decomposition is far more effective as a preconditioner than as a stand-alone solver, and until recently, only linear DD preconditioners have been available. An attractive framework to overcome this bottleneck and thus facilitate nonlinear preconditioning has emerged as the Additive Schwarz preconditioned inexact Newton (ASPIN) methods introduced by Cai and Keyes [13]. ASPIN has been shown to work well for several computational fluid dynamics problems, in that it reduces the computational cost, and can handle higher nonlinearities than standard Newton methods. Until now, no applications of ASPIN to porous media problems have been reported. Nonlinear preconditioning has however been considered in the context of mortar finite element methods [14], [15], [16], [17].

In this paper, we apply linear and nonlinear domain decomposition both as stand-alone solvers and as preconditioners for flow and transport in porous media. We consider two-phase flow problems with a range of difficulty that spans from stable displacement in homogeneous media to highly unstable displacement in channelized systems. To highlight nonlinear effects, the time stepping is fully implicit. The simulation results consistently show that the performance of ASPIN is better than that of linear preconditioning, by as much as 50–75% for the most challenging cases, depending on the number of subdomains used. Crucially, the performance of ASPIN is more robust than linear preconditioning in the face of increasing nonlinearities, indicating that ASPIN is the domain decomposition method of choice for challenging problems. Our main conclusion is thus that contrary to common practice, domain decomposition should be considered a preconditioner for the nonlinear system, not the linear system, for optimal performance and robustness.

The rest of this paper is organized as follows. In Section 2, we describe the model problem and discretization techniques used in this study. Section 3 gives an overview of domain decomposition methods and the specific methods considered here. In Section 4 we present numerical results from a comparative study of the different methods. The paper is concluded in Section 5.