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 socalled 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 nonnegative 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 largescale 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 pressuredependent 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. Pressuredependent 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 pressuredependent 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 DarcyForchheimer 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 nonnegative 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 singlefield 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 diffusiontype 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 metricbased 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 metricbased 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.
Optimizationbased 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 loworder finite elements, whose shape functions are nonnegative 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], advectiondiffusion 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 infinitedimensional variational inequalities. Subsequently the field of finitedimensional 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, finitedimensional variational inequalities have been utilized to enforce DMP and the nonnegativity of concentrations under advectiondiffusion 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 advectiondiffusion 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 pressuredependent viscosity, and possess the following attractive features:

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

The devised computational framework is equipped to handle nonlinear formulations and problems with nonselfadjoint operators.

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

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

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 NonLinear 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(2.1) 
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 secondorder 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
(2.2) 
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 twoterm 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:
(2.3) 
For convenience, we introduce the drag coefficient, which is defined as follows:
(2.4) 
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 pressuredependent viscosity can be written as follows:
(2.5a)  
(2.5b)  
(2.5c)  
(2.5d) 
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:
(2.6a)  
(2.6b) 
Note that equation (2.6a) is obtained by combining equations (2.5a) and (2.5b). Equation (2.6a) is a special case of a secondorder quasilinear elliptic partial differential equation.
A general secondorder quasilinear elliptic operator takes the following form:
(2.7) 
We define the coefficient matrix as follows:
(2.8) 
which, in indicial notation, takes the following form:
(2.9) 
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:

and are continuously differentiable with respect to and ,

is elliptic in ,

is nonincreasing with respect to for fixed , and

then we have in .
Proof.
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 nonnegative 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
(2.10) 
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:
(2.11) 
with the following properties:

and are continuously differentiable with respect to and ,

is elliptic in ,

is nonincreasing with respect to , and

.
Proof.
It is a straightforward computation to show that with the following choices:
(2.12a)  
(2.12b) 
equation (2.6a) can be written as
(2.13) 
implies that conditions (iii), (iv) and the second part of (i) are trivially satisfied. Using equation (2.2) we have
(2.14) 
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:
(2.15) 
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.
(Nonnegative pressures under the porous media model) Let be sufficiently smooth and . If the prescribed volumetric source is nonnegative in and the prescribed pressure on the boundary is nonnegative then the pressure in the entire domain is nonnegative. That is, if in and on then in .
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,
(2.20) 
Proof.
To show the right inequality, we take
(2.21) 
Since on , we have on . Moreover, the above choices for and imply that
(2.22) 
Lemma 2.2 and the comparison principle given by Theorem 2.1 imply that
(2.23) 
This further implies that
(2.24) 
One can thus conclude that .
To show the left inequality, we take
(2.25) 
which imply that
(2.26)  
(2.27) 
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 nonnegative 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 wellestablished 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 LadyzhenskayaBabǔskaBrezzi 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 lowestorder RaviartThomas space
[Raviart and Thomas, 1977]. It is wellknown 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 RaviartThomas 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 Multiscale 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:
(3.1a)  
(3.1b)  
(3.1c)  
(3.1d) 
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.2) 
3.1. Lowestorder RaviartThomas space
Given a simplex , the local RaviartThomas space of order is defined as follows [Raviart and Thomas, 1977; Bergamaschi et al., 1994]:
(3.3) 
where is the space of polynomials of degree and , as mentioned before, is the number of spatial dimensions. It is wellknown that the interpolations for the velocity and pressure fields under the RaviartThomas spaces of all orders satisfy the LBB infsup stability condition and thereby provide stable numerical solutions under the Galerkin weak formulation [Brezzi and Fortin, 1991].
In this study we employ the lowestorder RaviartThomas space for interpolation of velocity and pressure fields, which is the simplest and the most popular space among the class of RaviartThomas spaces. Under the lowestorder RaviartThomas 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,
(3.4) 
The finite dimensional subspaces and under for a triangle are defined as follows:
(3.5a)  
(3.5b) 
where is a triangulation on . These subspaces for tetrahedra are defined as follows:
(3.6a)  
(3.6b) 
where , in this case, is a tetrahedralization (i.e., 3D triangulation) on .
3.2. Variational Multiscale formulation
Variational Multiscale (VMS) is a computational paradigm to achieve enhanced stability of a given weak formulation [Hughes, 1995]. For a mixed formulation, say the Galerkin formulation, residualbased adjointtype 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
(3.7) 
In all our numerical simulations, we employ equalorder linear nodalbased interpolations for the pressure and velocity fields.
3.3. Nonlinear formulations
The pressure dependence of viscosity in the weak formulations turns the problem into a nonlinear problem. To solve such problems, we introduce the canonical form: Find and such that we have
(3.8) 
where is the residual expressed in semilinear form; the arguments to the left and right of the semicolon are nonlinear and linear, respectively. The semilinear form for the RT0 formulation takes the following form:
(3.9) 
The semilinear form for the VMS formulation can be written as follows:
(3.10) 
Newton’s method is employed to solve the nonlinear variational forms. Let the superscript () denote the current Newton or nonlinear 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:
(3.11) 
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:
(3.12) 
Likewise, the Jacobian for the VMS formulation reads:
(3.13) 
In each Newton iteration, we thus solve the following linear variational problem: Find and such that we have
(3.14) 
We obtain the solution in an iterative fashion using the following update equation until the residual meets the prescribed tolerance:
(3.15a)  
(3.15b) 
3.4. VI formulation in continuous setting
In order to enforce the bound constraints due to maximum principles and the nonnegative constraint we pose the problem as a variational inequality. To this end, we define the feasible solution space to be as follow:
(3.16) 
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
(3.17) 
4. Proposed Computational Framework
4.1. Solver methodology
Our proposed computational framework based on the two nonlinear finite element variational formulations result in saddlepoint problems, which are notoriously difficult to solve in a largescale 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 saddlepoint system.
The residual vector for the RT0 formulation can be written as:
(4.1a)  
(4.1b) 
where the subscripts and denote the velocity and pressure components respectively. Likewise, the residual vector for the VMS formulation is written as follows:
(4.2a)  
(4.2b) 
The components of the Jacobian matrices for equations (4.1) and (4.2), respectively, can be subdivided as follows:
(4.3a)  
(4.3b)  
(4.3c)  
(4.3d) 
and
(4.4a)  
(4.4b)  
(4.4c)  
(4.4d) 
Conceptually, the problem at hand is a 22 block matrix:
(4.5) 
which admits a full factorization of
(4.6) 
where
is the identity matrix and
(4.7) 
is the Schur complement. The inverse can therefore be written as
(4.8) 
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 masslumping of
to estimate
. That is,(4.9) 
to precondition the inner solver inverting . For this block we employ the multigrid Vcycle 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 Vcycle 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 degreesoffreedom by “
”. The componentwise inequalities are denoted by and . That is,(4.10a)  
(4.10b) 
The standard innerproduct in Euclidean spaces is denoted by . That is,
(4.11) 
Let and denote the discrete vector of unknowns for velocity and pressure respectively. The vector of all degreesoffreedom, denoted by , can be defined as:
(4.12) 
For convenience, let us also define the following functional as
(4.13) 
The VI formulation in the discrete setting is posed as a Mixed Complementarity Problem (MCP) [Kinderlehrer and Stampacchia, 2000]: Find such that for each
(4.14a)  
(4.14b)  
(4.14c) 
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 nonlinear complementarity problem, which is a special case of MCP. For details on nonlinear 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.15) 
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 pythonbased library that provides an automated system for the solution of partial differential equations using the finite element method. The MPIbased PETSc library is utilized as the parallel linear algebra backend. These solvers have been demonstrated to show good parallel scalability for largescale optimizationbased problems [Chang et al., 2017a].
The PETSc library [Balay et al., 2014] provides a wide array of solvers for finitedimensional VI’s. For example, two popular algorithms are the semismooth Newton (SS) [Luca et al., 1996; Munson et al., 2001] and Reducedspace activeset (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 advectiondiffusion 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:

Read in mesh, boundary conditions, and material properties.

Solve for with no constraints (call it ).

CONDITIONAL: If violates the discrete maximum principles:

Set as initial guess for SS method.

Solve for using SS method.

In the next section, all 2D problems will be conducted in serial on an Intel Xeon E52609v3 (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
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:
(5.1a)  
(5.1b) 
Through the method of manufactured solutions, by substituting equation (5.1) into equation (2.5) we obtain the following expression for the body force:
(5.2) 
where is given by equation (2.3). The boundary conditions for this problem are:
Parameter  Values 

1  
0 and 1  
1  
(5.3) 
Figure 1 provides a pictorial description of the problem, and Table LABEL:table:_hconv lists the parameters employed in the numerical simulation.
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 pressuredependent 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
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:
(5.4) 
where m denotes the base permeability and is a userdefined value that controls the level of anisotropy. Three values of are considered as shown in Table LABEL:table:_eps_cases.
Case ID  Value of 

1  
2  
3 
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:
(5.5) 
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 advectiondiffusion 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.53E001  36  2  6.26E001  1.58E+000  0.02 
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 
5.3. Circular reservoir
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 degreesoffreedom 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  VDOF  PDOF  TotalDOF  VDOF  PDOF  TotalDOF 
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 
Case ID  Units  

1  
2  
3 
Let Pas and the anisotropic permeability tensor be given by:
(5.6) 
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 
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.
Comments
There are no comments yet.