A scalable variational inequality approach for flow through porous media models with pressure-dependent viscosity

by   N. K. Mapakshi, et al.
University of Houston

Mathematical models for flow through porous media typically enjoy the so-called maximum principles, which place bounds on the pressure field. It is highly desirable to preserve these bounds on the pressure field in predictive numerical simulations, that is, one needs to satisfy discrete maximum principles (DMP). Unfortunately, many of the existing formulations for flow through porous media models do not satisfy DMP. This paper presents a robust, scalable numerical formulation based on variational inequalities (VI), to model non-linear flows through heterogeneous, anisotropic porous media without violating DMP. VI is an optimization technique that places bounds on the numerical solutions of partial differential equations. To crystallize the ideas, a modification to Darcy equations by taking into account pressure-dependent viscosity will be discretized using the lowest-order Raviart-Thomas (RT0) and Variational Multi-scale (VMS) finite element formulations. It will be shown that these formulations violate DMP, and, in fact, these violations increase with an increase in anisotropy. It will be shown that the proposed VI-based formulation provides a viable route to enforce DMP. Moreover, it will be shown that the proposed formulation is scalable, and can work with any numerical discretization and weak form. Parallel scalability on modern computational platforms will be illustrated through strong-scaling studies, which will prove the efficiency of the proposed formulation in a parallel setting. Algorithmic scalability as the problem size is scaled up will be demonstrated through novel static-scaling studies. The performed static-scaling studies can serve as a guide for users to be able to select an appropriate discretization for a given problem size.



There are no comments yet.


page 1

page 19

page 20

page 22

page 23

page 28

page 29

page 33


Numerical approximation of singular-degenerate parabolic stochastic PDEs

We study a general class of singular degenerate parabolic stochastic par...

Two mixed finite element formulations for the weak imposition of the Neumann boundary conditions for the Darcy flow

We propose two different discrete formulations for the weak imposition o...

A Variational Finite Element Discretization of Compressible Flow

We present a finite element variational integrator for compressible flow...

Free energy diminishing discretization of Darcy-Forchheimer flow in poroelastic media

In this paper, we develop a discretization for the non-linear coupled mo...

Adaptive and Pressure-Robust Discretization of Incompressible Pressure-Driven Phase-Field Fracture

In this work, we consider pressurized phase-field fracture problems in n...

Free surface flow through rigid porous media – An overview and comparison of formulations

In many applications free surface flow through rigid porous media has to...

Variational multiscale modeling with discretely divergence-free subscales

We introduce a residual-based stabilized formulation for incompressible ...
This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1. Introduction

The success of many current and emerging technological endeavors critically depend on a firm understanding and on the ability to control flows in heterogeneous, anisotropic porous media. These endeavors include geological carbon sequestration, geothermal systems, oil recovery, water purification systems, extraction of gas hydrates from tight shale; just to name a few. Modeling and predictive simulations play an important role in all these endeavors, and one has to overcome many numerical challenges to obtain accurate numerical solutions. It is beyond the scope of this paper to address all the major issues associated with the flow of fluids through porous media. Herein, we however address one of the main numerical challenges that is encountered in numerical modeling of flow through porous media with relevance to the mentioned applications.

Flow through porous media models typically enjoy the so-called maximum principles, which place bounds on the pressure field. These bounds depend on the prescribed data, which include boundary conditions, anisotropy of the porous media, body force, volumetric source, topology of the domain, and the regularity of the boundary. The non-negative constraint on the pressure (which basically implies the physical condition that a fluid subject to a flow in a porous medium cannot sustain a “suction” by itself) can be shown to be a special case of the classical maximum principle. It is imperative that these bounds on the pressure field are preserved in a predictive numerical simulation; that is, one needs to satisfy maximum principles in the discrete setting. The discrete version of maximum principles is commonly referred to as discrete maximum principles (DMP). It becomes even more crucial for those flow models in which the material properties depend on the pressure; for example, the case in which the viscosity of the fluid depends on the pressure in the fluid, as a violation of DMP can amplify errors in the solution fields. Unfortunately, many of the commonly used mixed finite element formulations for flow through porous media models do not satisfy DMP, which will be shown in the subsequent sections. Moreover, the problems pertaining to flow through porous media, especially the ones encountered in subsurface modeling, are highly nonlinear and large-scale in nature. Thus, one needs to develop numerical formulations that are scalable in an algorithmic and parallel sense in addition to satisfying DMP.

This paper presents a new, scalable numerical formulation based on variational inequalities (VI) that enforces discrete maximum principles for nonlinear flow through porous media models by taking into account heterogeneity, anisotropic permeability and pressure-dependent viscosity.

1.1. A review of related prior works

In order to bring out clearly the contributions made in this paper and the approach taken by us, we provide a brief discussion on prior works with respect to three aspects.

1.1.1. Pressure-dependent viscosity.

The classical Darcy model [Darcy, 1856], which is the most popular flow through porous media model, assumes the viscosity of the fluid to be a constant, and in particular, the model assumes that the coefficient of viscosity is independent of the pressure in the fluid [Nakshatrala and Rajagopal, 2011]. But there is abundant experimental evidence that the viscosity of liquids, especially organic liquids, depends on the pressure [Bridgman, 1931]. More importantly, the dependence of viscosity on pressure for organic liquids is exponential [Barus, 1893]. Since then several studies have developed mathematical models that take into account the dependence of viscosity on pressure, and established the existence of solutions for the resulting governing equations [Málek et al., 2002; Hron et al., 2003; Franta et al., 2005; Bulíček et al., 2007]. A work that is relevant to this paper is by [Nakshatrala and Rajagopal, 2011] who derived a modification to the Darcy model using the mixture theory by taking into account the pressure-dependent viscosity. They have also developed a stabilized formulation for the resulting equations using the variational multiscale paradigm [Hughes, 1995], and have shown, using numerical simulations, that the dependence of viscosity on pressure has a significant effect on both qualitative and quantitative nature of the solution fields. Later, [Nakshatrala and Turner, 2013] have presented a stabilized mixed formulation based on Picard linearization and laid down the differences in the predictions for enhanced oil recovery and carbon sequestration when the pressure dependence on viscosity is considered against when not considered. Recently, [Chang et al., 2017c] have extended the pressure dependence to the Darcy-Forchheimer model, and demonstrated how the dependence of the drag coefficient on pressure differs significantly from when it depends on velocity. However, all of these studies considered isotropic permeability, and did not address the violations of maximum principles and the non-negative constraint on the pressure field.

1.1.2. Anisotropy, violations of DMP and numerical techniques to enforce DMP

[Varga, 1966] was the first to address DMP, and the study was restricted to the finite difference method applied on the Poisson’s equation. [Ciarlet and Raviart, 1973] were the first to address DMP in the context of the finite element method. Their study revealed that the single-field Galerkin formulation for solving the Poisson’s equation, in general, does not satisfy DMP. They also obtained sufficient conditions which are in the form of restrictions on the mesh (e.g., all the angles of a triangular element to be acute) to meet DMP. Subsequent studies have found that these mesh restrictions, which have been derived for isotropic diffusion equations, are not sufficient when one considers anisotropic diffusion equations or other processes like advection and reactions; for example, see [Mudunuru and Nakshatrala, 2017] and references therein.

In the last decade, several approaches have been developed to enforce DMP on general computational grids for anisotropic diffusion-type equations under the finite element method. Some of the notable approaches are based on either constrained optimization techniques [Liska and Shashkov, 2008; Nagarajan and Nakshatrala, 2011; Mudunuru et al., 2015], placing anisotropic metric-based restrictions on the mesh [Huang and Wang, 2015; Mudunuru and Nakshatrala, 2017], or altering the formulations at the continuum setting [Pal et al., 2016]. Placing restrictions on meshes is not a viable approach for applications involving flow through porous media, as the computational domains are complex and it is not practical or even possible to generate metric-based meshes that satisfy DMP. The approach of altering formulations at the continuum setting is not a viable route either for porous media applications, as one has to deal with a hierarchy of models with multiple constituents in such applications and there is no trivial way of altering formulations at the continuum level to meet DMP.

Optimization-based methods based on quadratic programming have been successfully employed to develop formulations for anisotropic diffusion equations, for example, see [Nakshatrala and Valocchi, 2009]. The key behind these methods has been to construct an objective function in quadratic form; employ low-order finite elements, whose shape functions are non-negative within each element; and enforce the constraints arising from DMP as explicit bound constraints on the nodal quantities. This approach has also been extended to transient problems [Nakshatrala et al., 2016], advection-diffusion equation [Mudunuru and Nakshatrala, 2016], diffusion with fast reactions [Nakshatrala et al., 2013], and parallel environments [Chang et al., 2017a]. However, it needs to be emphasized that this approach, which is based on quadratic programming, requires the bilinear form in the weak form to be symmetric, which is not the case with many porous media models and weak formulations. More importantly, all the mentioned studies considered linear equations arising from transport problems.

1.1.3. Variational inequality based techniques.

Variational inequalities arise quite naturally in various branches of mechanics [Kikuchi and Oden, 1988; Hlavacek et al., 2012; Rodrigues, 1987; Han and Reddy, 2012]. In fact, the whole field of variational inequalities grew from a problem in mechanics which was posed by [Signorini, 1933]. This problem, which is popularly referred to as the Signorini problem, is about finding static equilibrium configurations of a linear elastic body resting on a rigid smooth surface [Signorini, 1959]. A work that is more directly related to this paper is by [Chipot, 2012], who has employed variational inequalities to address some class of problems that arise in studies on flow through porous media, specifically the dam problem, and established mathematical properties like existence and uniqueness of solutions. The Signorini problem and the treatment of the dam problem are examples of infinite-dimensional variational inequalities. Subsequently the field of finite-dimensional variational inequalities has been developed [Facchinei and Pang, 2003], and this field has eventually found its way into the mainstream numerical optimization [Kinderlehrer and Stampacchia, 2000; Ulbrich, 2011]. Recently, finite-dimensional variational inequalities have been utilized to enforce DMP and the non-negativity of concentrations under advection-diffusion equations [Chang and Nakshatrala, 2017]. This formulation does not require the bilinear form of the underlying weak formulation to be symmetric. However, it needs to be emphasized that advection-diffusion equations, which arise in transport problems, are linear and are in terms of a single unknown field. In this paper, we address flow problems, and the governing equations are in terms of two fields (i.e., the velocity and pressure) and are nonlinear.

1.2. Our approach, salient features and an outline of the paper

We develop a numerical formulation based on variational inequalities for modeling flow through porous media that accounts for anisotropy and pressure-dependent viscosity, and possess the following attractive features:

  1. The proposed framework can handle any numerical discretization and weak formulation.

  2. The devised computational framework is equipped to handle non-linear formulations and problems with non-self-adjoint operators.

  3. Maximum principles are satisfied under this framework even when anisotropy is present.

  4. A computer implementation of the proposed framework can seamlessly leverage on the state-of-the-art software and algorithms that are currently available for high performance computing.

  5. The implementation outlined in this paper has excellent scalability both in the algorithmic and parallel sense.

All the aforementioned features of the proposed formulation will be illustrated in the subsequent sections. The rest of this paper is organized as follows. In Section 2 we present the governing equations for the modified Darcy flow and describe discrete maximum principles. In Section 3 we provide various mixed and nonlinear formulations used for this study. In Section 4, we lay down the solver methodology and outline of the computer implementation for the proposed VI based framework. In Section 5, we present numerical results illustrating effectiveness and scalability of the proposed framework. Concluding remarks are made in Section 6.

2. A Non-Linear Model for Flow Through Porous Media

Consider a porous domain denoted by , where “” denotes the number of spatial dimensions. The boundary of the domain will be denoted by , where an overline denotes the set closure. A spatial point will be denoted by . The divergence and gradient operators with respect to are, respectively, denoted by and . The unit outward normal to the boundary is denoted by

. The discharge velocity vector field and the pressure scalar field are denoted by

and , respectively. The boundary is divided into two parts:  and , such that we have


denotes that part of the boundary on which the normal component of the velocity is prescribed. is that part of the boundary on which pressure is prescribed. The permeability of the porous medium will be denoted by

, which is a second-order tensor. It is assumed that the permeability tensor is positive definite and symmetric. The density of the fluid is denoted by

. The coefficient of (dynamic) viscosity of the fluid is denoted by .

As mentioned earlier, for most organic liquids, the dependence of pressure on viscosity is exponential. That is, mathematically we have


where is the Barus coefficient, which has to be obtained experimentally and the values of this coefficient for various liquids can be found in [Bridgman, 1931]. Since the pressures that we deal in this paper are relatively small but not sufficiently small enough to neglect the pressure dependence of viscosity, we take a two-term Taylor expansion of the Barus formula (given by equation (2.2)). The two term Taylor expansion, which will be employed in all the numerical experiments in this paper, takes the following mathematical form:


For convenience, we introduce the drag coefficient, which is defined as follows:


Since the viscosity depends on the pressure and the permeability explicitly depends on the spatial coordinates, the drag coefficient will explicitly depend on both the pressure and .

The governing equations for flow through porous media by taking into account the pressure-dependent viscosity can be written as follows:


where denotes the prescribed normal component of velocity on the boundary, denotes the prescribed pressure on the boundary, is the specific body force, and is the prescribed volumetric source, all of which are functions of . It should be noted that one can recover the classical Darcy equations by setting , which makes the drag coefficient to be independent of the pressure. A systematic derivation of the above governing equations under the theory of interacting continua can be found in [Nakshatrala and Rajagopal, 2011].

2.1. A mathematical interlude

For a mathematical treatment of the abstract boundary value problem (2.5), we assume to have pressure boundary conditions on the entire boundary (i.e., ). We also rewrite the governing equations solely in terms of the pressure as follows:


Note that equation (2.6a) is obtained by combining equations (2.5a) and (2.5b). Equation (2.6a) is a special case of a second-order quasi-linear elliptic partial differential equation.

A general second-order quasi-linear elliptic operator takes the following form:


We define the coefficient matrix as follows:


which, in indicial notation, takes the following form:


Note that the entries of the coefficient matrix need not be constants. The operator is said to be elliptic in if the coefficient matrix is positive definite for all , and . From the theory of partial differential equations, this operator is known to satisfy the following important property:

Theorem 2.1.

(Comparison principle in a general setting) Let which satisfy and in and on . If the following conditions are met:

  1. and are continuously differentiable with respect to and ,

  2. is elliptic in ,

  3. is non-increasing with respect to for fixed , and

then we have in .


A mathematical proof can be found in [Gilbarg and Trudinger, 2001, Theorem 10.7]. ∎

We now show that the solutions of the porous media model satisfy the non-negative constraint and the maximum principle. To this end, we assume that the body force is conservative. That is, there exists a scalar field such that

Lemma 2.2.

The porous media model, given by equations (2.6a), along with the pressure dependent viscosity given by equation (2.2), can be put into the following form:


with the following properties:

  1. and are continuously differentiable with respect to and ,

  2. is elliptic in ,

  3. is non-increasing with respect to , and

  4. .


It is a straightforward computation to show that with the following choices:


equation (2.6a) can be written as


implies that conditions (iii), (iv) and the second part of (i) are trivially satisfied. Using equation (2.2) we have


Clearly, is continuously differentiable with respect to and , which implies that the first part of condition (i) is satisfied. The coefficient matrix can be written as follows:


The positive definiteness of the permeability tensor, , and imply that the coefficient matrix is positive definite for all , and ; which further implies that condition (ii) is met. ∎

Theorem 2.3.

(Non-negative pressures under the porous media model) Let be sufficiently smooth and . If the prescribed volumetric source is non-negative in and the prescribed pressure on the boundary is non-negative then the pressure in the entire domain is non-negative. That is, if in and on then in .


Choose and . We then have


If on we have


Using Lemma 2.2 and the comparison principle given by Theorem 2.1, we conclude that


This further implies that


which completes the proof. ∎

Theorem 2.4.

(Maximum principle for the porous media model) If the prescribed volumetric source is zero (i.e., ) in then the maximum and minimum pressures occur on the boundary. That is,


To show the right inequality, we take


Since on , we have on . Moreover, the above choices for and imply that


Lemma 2.2 and the comparison principle given by Theorem 2.1 imply that


This further implies that


One can thus conclude that .

To show the left inequality, we take


which imply that


By again appealing to the comparison principle, we conclude that in . ∎

The existing numerical discretizations for flow through porous media models do not produce solutions that satisfy the aforementioned mathematical properties for anisotropic porous domains. Thus, the central aim of this paper is to develop a computational framework for nonlinear models for flow through porous media that satisfies the maximum principle and ensures non-negative solutions for the pressure. This will be achieved by combining mixed finite element methods and variational inequalities.

3. Mixed Formulations

In our study, we employ two well-established finite element formulations which achieve discrete stability differently. The stability of a mixed formulation in the discrete setting will be primarily dictated by the famous Ladyzhenskaya-Babǔska-Brezzi stability condition [Babuška, 1973; Brezzi and Fortin, 1991]

. The first formulation is the classical mixed formulation (also known as the Galerkin weak formulation) but the interpolations for the velocity and pressure fields are based on the lowest-order Raviart-Thomas space

[Raviart and Thomas, 1977]. It is well-known that an arbitrary combination of interpolation functions for the velocity and pressure fields under the classical mixed formulation need not satisfy the LBB condition, and hence may not be stable [Brezzi and Fortin, 1991; Brezzi et al., 2008]. The Raviart-Thomas spaces place restrictions on the interpolations for the velocity and pressures fields to satisfy the LBB condition, and thus making the classical mixed formulation stable. The second formulation is the Variational Multi-scale formulation [Nakshatrala and Rajagopal, 2011], which is a stabilized mixed formulation that augments the Galerkin weak formulation with stabilization terms to circumvent the LBB condition. A nice discussion on the two classes of mixed formulation, which differ in the way they handle the LBB condition (i.e., satisfying vs. circumventing), can be found in [Franca and Hughes, 1988].

The weak forms of the aforementioned two formulations will form the basis for the associated variational inequalities. To this end, the following function spaces will be employed in the rest of the paper:


where is a standard Sobolov space [Brezzi and Fortin, 1991] and is set of all square integrable functions on . The Galerkin weak formulation for the governing equations (2.5) reads: Find and such that we have


3.1. Lowest-order Raviart-Thomas space

Given a simplex , the local Raviart-Thomas space of order is defined as follows [Raviart and Thomas, 1977; Bergamaschi et al., 1994]:


where is the space of polynomials of degree and , as mentioned before, is the number of spatial dimensions. It is well-known that the interpolations for the velocity and pressure fields under the Raviart-Thomas spaces of all orders satisfy the LBB inf-sup stability condition and thereby provide stable numerical solutions under the Galerkin weak formulation [Brezzi and Fortin, 1991].

In this study we employ the lowest-order Raviart-Thomas space for interpolation of velocity and pressure fields, which is the simplest and the most popular space among the class of Raviart-Thomas spaces. Under the lowest-order Raviart-Thomas space, the pressure is constant within an element and the fluxes are evaluated at the midpoint of each edge in 2D or at the barycenter of each face in 3D. Mathematically,


The finite dimensional subspaces and under for a triangle are defined as follows:


where is a triangulation on . These subspaces for tetrahedra are defined as follows:


where , in this case, is a tetrahedralization (i.e., 3D triangulation) on .

3.2. Variational Multi-scale formulation

Variational Multi-scale (VMS) is a computational paradigm to achieve enhanced stability of a given weak formulation [Hughes, 1995]. For a mixed formulation, say the Galerkin formulation, residual-based adjoint-type stabilization terms are added to circumvent the LBB condition and achieve stability. [Nakshatrala and Rajagopal, 2011] have successfully employed the VMS paradigm to develop a stabilized mixed formulation for the isotropic version of the porous media model outlined in Section 2. The weak form under the VMS formulation for governing equations (2.5a)–(2.5d) reads: Find and such that we have


In all our numerical simulations, we employ equal-order linear nodal-based interpolations for the pressure and velocity fields.

3.3. Non-linear formulations

The pressure dependence of viscosity in the weak formulations turns the problem into a non-linear problem. To solve such problems, we introduce the canonical form: Find and such that we have


where is the residual expressed in semi-linear form; the arguments to the left and right of the semicolon are non-linear and linear, respectively. The semi-linear form for the RT0 formulation takes the following form:


The semi-linear form for the VMS formulation can be written as follows:


Newton’s method is employed to solve the non-linear variational forms. Let the superscript () denote the current Newton or non-linear iteration. The Jacobian is computed by taking the Gteaux variation of the residual at and in the directions of and respectively. Formally, this is derived by computing:


provided the limit exists. For further details on the Gâteaux variation see [Spivak, 1997; Holzapfel, 2000; Glowinski, 2008]. Following through with the calculation above yields the following Jacobian under the RT0 formulation:


Likewise, the Jacobian for the VMS formulation reads:


In each Newton iteration, we thus solve the following linear variational problem: Find and such that we have


We obtain the solution in an iterative fashion using the following update equation until the residual meets the prescribed tolerance:


3.4. VI formulation in continuous setting

In order to enforce the bound constraints due to maximum principles and the non-negative constraint we pose the problem as a variational inequality. To this end, we define the feasible solution space to be as follow:


In a specific problem, if there is no restriction on the lower bound of the pressure then one can set . Similarly, one can set if there is no upper bound on the pressure. The proposed variational inequality in the continuous setting reads: Find such that we have


4. Proposed Computational Framework

4.1. Solver methodology

Our proposed computational framework based on the two non-linear finite element variational formulations result in saddle-point problems, which are notoriously difficult to solve in a large-scale setting. Several classes of iterative solvers and preconditioning strategies exist for these types of problems [Benzi et al., 2005; Elman et al., 2006; Murphy et al., 2000]. One could alternatively employ hybridization techniques [Cockburn et al., 2009] which introduces Lagrange multipliers which can also significantly reduce the difficulty of solving such problems. However, in this study, we employ a Schur complement approach to precondition the saddle-point system.

The residual vector for the RT0 formulation can be written as:


where the subscripts and denote the velocity and pressure components respectively. Likewise, the residual vector for the VMS formulation is written as follows:


The components of the Jacobian matrices for equations (4.1) and (4.2), respectively, can be subdivided as follows:




Conceptually, the problem at hand is a 22 block matrix:


which admits a full factorization of



is the identity matrix and


is the Schur complement. The inverse can therefore be written as


The task at hand is to approximate and . Since is a mass matrix for the Darcy equation, we can invert it using the ILU(0) (incomplete lower upper) solver. We employ a diagonal mass-lumping of

to estimate

. That is,


to precondition the inner solver inverting . For this block we employ the multi-grid V-cycle on from the HYPRE BoomerAMG package ([Falgout, 2006]). As discussed in Appendix B of [Chang and Nakshatrala, 2017], the and , only a single sweep of ILU(0) and HYPRE’s V-cycle is needed for the and matrices, respectively, and the GMRES method is employed to solve the entire block system.

4.2. Variational inequality approach

We denote the total number of degrees-of-freedom by “

”. The component-wise inequalities are denoted by and . That is,


The standard inner-product in Euclidean spaces is denoted by . That is,


Let and denote the discrete vector of unknowns for velocity and pressure respectively. The vector of all degrees-of-freedom, denoted by , can be defined as:


For convenience, let us also define the following functional as


The VI formulation in the discrete setting is posed as a Mixed Complementarity Problem (MCP) [Kinderlehrer and Stampacchia, 2000]: Find such that for each


where and , respectively, denote the minimum and maximum values for pressure and velocity. The constraints for pressure ( and ) are provided by the maximum principle, whereas the minimum and maximum constraints for each directional component of velocity are and respectively.

If one has only lower bound constraints due to the presence of a positive pressure source (i.e., , , and ), then the VI reduces to a non-linear complementarity problem, which is a special case of MCP. For details on non-linear complementarity problems, see [Facchinei and Pang, 2003]. Note that the feasible region, which is restricted by the bound constraints, form a parallelepiped, which is a convex set [Boyd and Vandenberghe, 2004].

Let the feasible region be a convex subset of . In our case, the feasible region is restricted by constraints which are in the form of finite number of linear equalities and inequalities. This makes the feasible region to be a polyhedron, which is a convex set [Boyd and Vandenberghe, 2004]. It should be noted that bound constraints are a special case of linear inequalities. With this machinery at our disposal, one can pose the second formulation based on variational inequalities, which reads: Find such that we have


4.3. Computer implementation

In this paper, we implement the proposed variational inequality based computational framework using the Firedrake Project [Rathgeber et al., 2016; Luporini et al., 2016, 2015]. It is a python-based library that provides an automated system for the solution of partial differential equations using the finite element method. The MPI-based PETSc library is utilized as the parallel linear algebra back-end. These solvers have been demonstrated to show good parallel scalability for large-scale optimization-based problems [Chang et al., 2017a].

The PETSc library [Balay et al., 2014] provides a wide array of solvers for finite-dimensional VI’s. For example, two popular algorithms are the semi-smooth Newton (SS) [Luca et al., 1996; Munson et al., 2001] and Reduced-space active-set (RS) [Benson and Munson, 2006] methods. It has been shown in [Benson and Munson, 2006] that the performance of SS and RS methods are application dependent and in [Chang and Nakshatrala, 2017] that the RS method demonstrates better solver convergence for advection-diffusion type equations. However, preliminary results (not shown in the paper) suggest that it is in fact the SS method that performs better for the nonlinear flow model. Thus, we propose the following algorithm:

  1. Read in mesh, boundary conditions, and material properties.

  2. Solve for with no constraints (call it ).

  3. CONDITIONAL: If violates the discrete maximum principles:

    1. Set as initial guess for SS method.

    2. Solve for using SS method.

In the next section, all 2D problems will be conducted in serial on an Intel Xeon E5-2609v3 (Haswell) processor, and the 3D problem will be conducted in parallel on an Intel Xeon Phi 7250 (Knights Landing) processor.

5. Representative Numerical Results

5.1. - convergence study

Figure 1. Problem statement for the -convergence study, showing the computational domain, pressure and velocity functions and the boundary conditions.

We first perform an -convergence study on the mixed formulations to verify that the Firedrake project library and proposed solver methodologies are converging schemes. A finite element solution is said to be converging if the difference between the exact and numerical solutions decreases with the mesh refinement. Consider a unit square to be the computational domain with the following expressions for the pressure and velocity fields:


Through the method of manufactured solutions, by substituting equation (5.1) into equation (2.5) we obtain the following expression for the body force:


where is given by equation (2.3). The boundary conditions for this problem are:

Parameter Values
0 and 1
Table 1. User defined parameters for -convergence study

Figure 1 provides a pictorial description of the problem, and Table LABEL:table:_h-conv lists the parameters employed in the numerical simulation.

Figure 2. Convergence plots of errors for RT0 and VMS formulations for the standard Darcy and modified Darcy models.

Figure 2 provides a comparison between error norm convergence rates for RT0 and VMS. Theoretical convergence rates of both velocity and pressure under RT0 is unity whereas the convergence rates for these two fields under VMS is two. We see that the slopes for the standard Darcy model (where ) are similar to the theoretical slopes, which verifies the convergence of the Firedrake project’s finite element framework. Furthermore, extending the mixed formulations to the modified Darcy model with pressure-dependent viscosity (where ) has similar convergence. The studies performed so far suggest that the Firedrake project is a suitable software package for conducting finite element simulations, and we now examine scenarios where VI is needed to enforce maximum principles.

5.2. Square reservoir

Figure 3. Pictorial description of square reservoir problem. The computational domain, boundary conditions, and volumetric source are shown.

This 2D heterogeneous problem aims to illustrate not only the effectiveness of the proposed VI based framework to ensure DMP for the pressure field but also how levels of heterogeneity and anisotropy affect the overall computational effort. Consider a square reservoir on a domain : = (0 m, 100 m) (0 m, 100 m) with the following anisotropic heterogeneous permeability tensor:


where m denotes the base permeability and is a user-defined value that controls the level of anisotropy. Three values of are considered as shown in Table LABEL:table:_eps_cases.

Case ID Value of
Table 2. Square reservoir: different values of used for this study.
(a) No VI ()
(b) With VI ()
(c) No VI ()
(d) With VI ()
(e) No VI ()
(f) With VI ()
Figure 4. Square reservoir: pressure contours for RT0 formulation before (left) and after (right) VI. The white spaces are representative of DMP violations.
(a) No VI ()
(b) With VI ()
(c) No VI ()
(d) With VI ()
(e) No VI ()
(f) With VI ()
Figure 5. Square reservoir: pressure contours for VMS formulation before (left) and after (right) VI. The white spaces are representative of DMP violations.

We assume that Pas and Pa. A constant pressure of 101325 Pa (1 atm) is applied on the entire boundary for both mixed formulations, see Figure 3 for a pictorial description. Let and the volumetric source be given by:


Since a positive forcing function is present, only the lower bound constraint is enforced in the VI framework. Both RT0 and VMS formulations are employed for the numerical discretization, with triangular elements of -size = 1 m. We perform this study for three cases of as listed in Table LABEL:table:_eps_cases.

Figures 4 and 5 show the pressure contours of the RT0 and VMS formulations, respectively, for different values of . The white regions are representative of DMP violations of pressure. Moreover, the effect of enforcing the VI framework over the RT0 and VMS formulations is also shown. These figures demonstrate that the VI framework is capable of enforcing the lower bound constraints for pressure values. It is also interesting to note that a smaller results in more violations regardless of the finite element discretization. Tables 3 and 4 illustrate the effect of on both the number of violating cells as well as solver performance. Small values of make the systems of equation much harder to solve as both time to solution and number of solver iterations increase. However, it can be seen that the additional computational cost from the VI solver is not significant; it was shown in [Chang and Nakshatrala, 2017] that the same computational framework with the RS method for advection-diffusion equations increased time to solution by up to a factor of 20 whereas for this particular problem the SS method increased time to solution by no more than 66 percent.

In a mixed formulation, altering the pressures or velocities may have a direct impact over their counterparts, but such changes to the velocity imposed by the VI framework do not have a negative impact on the overall numerical accuracy. The velocities under the VI framework for the RT0 and VMS formulations are shown in Figures 6 and 7, respectively. Absolute differences between velocity fields obtained from the mixed formulations without the imposed VI framework and the velocity fields obtained from the mixed formulations under the imposed VI framework are shown in those figures. These figures suggest that enforcement of the DMP on pressure does impact the velocities. Some formulations like RT0 have larger differences in velocity magnitudes whereas for others like VMS, it may not be as large. However, any such change to the velocity field can have a significant impact on subsurface transport especially for long periods of time.

Case RT0 VI over RT0 Total % Vio-
ID KSP SNES Time KSP SNES Time Time lations
1 2059 8 1.98E+001 664 3 4.53E+000 2.43E+001 50.54
2 227 4 2.63E+000 125 2 1.54E+000 4.18E+000 2.06
3 53 3 9.53E-001 36 2 6.26E-001 1.58E+000 0.02
Table 3. Square reservoir: computational results and initial violations for RT0.
Case VMS VI over VMS Total % Vio-
ID KSP SNES Time KSP SNES Time Time lations
1 608 5 1.75E+001 47 2 3.08E+000 2.06E+001 55.67
2 185 4 8.31E+000 71 2 3.80E+000 1.21E+001 2.37
3 51 3 5.00E+000 32 2 3.29E+000 8.29E+000 1.44
Table 4. Square reservoir: computational results and initial violations for VMS.
(a) VI velocity ()
(b) Absolute diff. ()
(c) VI velocity ()
(d) Absolute diff. ()
(e) VI velocity ()
(f) Absolute diff. ()
Figure 6. Square reservoir: velocity profiles for VI based framework imposed over RT0 formulation (left) and the absolute differences between the non-VI and VI velocities (right).
(a) VI velocity ()
(b) Absolute diff. ()
(c) VI velocity ()
(d) Absolute diff. ()
(e) VI velocity ()
(f) Absolute diff. ()
Figure 7. Square reservoir: velocity profiles for VI based framework imposed over VMS formulation (left) and the absolute differences between the non-VI and VI velocities (right).

5.3. Circular reservoir

Figure 8. Pictorial description of the circular reservoir problem, showing the computational domain and boundary conditions.

For this problem, the task is to study the algorithmic scalability of the proposed computational framework when different unstructured grids and different coefficients are employed. Consider the circular reservoir shown in Figure 8 which has a 100 m outer radius and an inner circular borehole radius of 1 m. Neither the specific body force nor the volumetric force is present for this problem so the VI framework imposes both lower bound and upper bound constraints. The lower bound or outer boundary is maintained at the atmospheric pressure ( Pa). At the injection hole, a constant pressure of 100 atm is maintained ( Pa) and serves as the upper bound constraint. The different meshes used for this study, as well as the corresponding numbers of degrees-of-freedom for the RT0 and VMS formulations, are provided in Table 5. Each mesh ID and mixed formulation shall be simulated with various values of provided in Table 6.

Mesh No.of No.of RT0 VMS
ID nodes elements V-DOF P-DOF Total-DOF V-DOF P-DOF Total-DOF
1 1379 2758 4061 2682 6743 2758 1379 4137
2 2485 4970 7744 5128 12872 5232 2616 7848
3 4822 9644 14326 9504 23830 9644 4822 14466
4 10356 20712 30858 20502 51360 20712 10356 31068
5 17673 35346 52741 35068 87809 35346 17673 53019
6 26926 53852 80440 53514 133954 53852 26926 80778
7 42911 85822 128301 85390 213691 85822 42911 128733
Table 5. Hierarchy of meshes for the circular reservoir problem.
Case ID Units
Table 6. Circular reservoir: values used for this study.

Let Pas and the anisotropic permeability tensor be given by:

Mesh ID % Violations: RT0 % Violations: VMS
Case 1 Case 2 Case 3 Case 1 Case 2 Case 3
1 56.52 57.34 62.60 48.77 49.03 49.89
2 53.33 54.43 59.89 46.04 46.20 47.11
3 54.63 55.66 59.04 45.01 47.92 48.98
4 53.19 53.92 56.73 47.95 50.60 51.55
5 51.06 51.60 54.46 51.50 46.29 50.13
6 52.85 51.87 53.01 47.25 48.81 52.58
7 52.00 50.02 51.65 47.52 36.12 49.62
Table 7. Circular reservoir: percentage of DMP violations for RT0 and VMS formulations.

where . Figures 9 and 10 show the pressures of RT0 and VMS formulations, respectively, for different values of before and after VI is imposed. It can be seen from these figures that DMP violations occur regardless of the used, and again the computational framework fixes these violations. Graphical representation of these violations are shown only for the first mesh (Mesh ID: 1), but detailed results concerning the violations for all other meshes are provided in Table 7.

(a) No VI (