Nonlinear nonlocal multicontinua upscaling framework and its applications
Abstract
In this paper, we discuss multiscale methods for nonlinear problems. We use recently developed multiscale concepts for linear problems and extend them to nonlinear problems. The main idea of these approaches is to use local constraints and solve problems in oversampled regions for constructing macroscopic equations. These techniques are intended for problems without scale separation and high contrast, which often occur in applications. For linear problems, the local solutions with constraints are used as basis functions. This technique is called Constraint Energy Minimizing Generalized Multiscale Finite Element Method (CEMGMsFEM). GMsFEM identifies macroscopic quantities based on rigorous analysis. In corresponding upscaling methods, the multiscale basis functions are selected such that the degrees of freedom have physical meanings, such as averages of the solution on each continuum.
This paper extends the linear concepts to nonlinear problems, where the local problems are nonlinear. The main concept consists of: (1) identifying macroscopic quantities; (2) constructing appropriate oversampled local problems with coarsegrid constraints; (3) formulating macroscopic equations. We consider two types of approaches. In the first approach, the solutions of local problems are used as basis functions (in a linear fashion) to solve nonlinear problems. This approach is simple to implement; however, it lacks the nonlinear interpolation, which we present in our second approach. In this approach, the local solutions are used as a nonlinear forward map from local averages (constraints) of the solution in oversampling region. This local finegrid solution is further used to formulate the coarsegrid problem. Both approaches are discussed on several examples and applied to singlephase and twophase flow problems, which are challenging because of convectiondominated nature of the concentration equation. The numerical results show that we can achieve good accuracy using our new concepts for these complex problems.
1 Introduction
As multiscale methods for linear equations are getting matured, their extension to difficult nonlinear problems remain challenging. Many nonlinear problems have multiscale nature due to spatial and temporal scales. For example, the dynamics of multiphase flow and transport in heterogeneous media varies over multiple space and time scales. In the past, many wellknown linear and nonlinear upscaling tools have been developed (e.g., [2, 27, 3, 26, 15, 8, 10, 23, 1, 21, 35, 36, 49, 47, 4, 43, 14, 17, 12, 52, 5, 13, 40, 48, 55]). Singlephase upscaling methods include permeability upscaling [22, 54, 9, 44] and many multiscale techniques [2, 26, 15, 49, 47, 4, 43]). Nonlinear upscaling methods, e.g., known as pseudorelative permeability approach [9, 45, 7], computes nonlinear relative permeability functions based on single cell twophase flow computations. It is known that these nonlinear approaches lack robustness and they are processes dependent [24, 25]. To overcome these difficulties, one needs a better understanding of nonlinear upscaling methods, which is our goal.
Nonlinear upscaling methods can be traced back to nonlinear homogenization [51, 30]. The main idea of nonlinear homogenization is to formulate coarsegrid equations based on nonlinear local problems formulated in each coarse block. Some examples include pLaplacian, pseudoelliptic equations, and parabolic equations [30]. In these approaches, the local problems are solved with periodic boundary conditions and some constraints on averages of the solutions of gradients. In nonperiodic cases, these methods are extended by solving local problems in coarse block subject to some boundary conditions. Because of nonlinearity, these local problems need to be solved for all possible average values, which can make the computations expensive. One can compute the upscaled fluxes onthefly using the values at previous iteration or previous time. These approaches are limited to problems without high contrast and scale separation. Our goal is to extend these approaches to problems with high contrast and nonseparable scales. The main novel components of our approach is introducing multiple macroscopic variables for each coarsegrid block, formulating appropriate local constrained problems that determine the downscaling map; formulating macroscopic equations.
In computational mechanics literature, there has been a great deal of research dedicated to nonlinear upscaling methods, which include generalized continuum theories (e.g., [32]), computational continua framework (e.g., [39]), and other approaches. Multiscale enrichment method based on the partition of unity ([41, 42]) combine the linear and nonlinear homogenization theory with the partition of unity method to handle the problems with inseparable fine and coarse scales. Computational continua ([39, 33]), which use nonlocal quadrature to couple the coarse scale system stated on a unions of some disjoint computational unit cells, are introduced for nonscaleseparation heterogeneous media. In [38, 37, 34], the authors further enhance the method by combining the computational continua with model reduction technique. Originally, the computational continua approaches were developed to overcome the limitations of higherorder homogenization models and generalized continuum theories, namely, the need for higherorder finite element continuity, additional degrees of freedom, and nonclassical boundary conditions. Though there are some similarities, our proposed approaches differ from computational continua framework. Our approach explores the localization of highcontrast (with respect to the coarsemesh size) features by introducing additional degrees of freedom (continua) and using oversampled regions. In future, we plan to use some ingredients of computational continua approaches in improving our approach and in making it more applicable to computational materials.
To design the new upscaled model, we use the concept of nonlocal multicontinuum (NLMC). The main idea of NLMC derives from Constraint Energy Minimizing Generalized Multiscale Finite Element Method (CEMGMsFEM) [19, 20, 16]. CEMGMsFEM identifies the degrees of freedom which can not be localized and then construct multiscale basis functions with support in the oversampled regions. However, the degrees of freedom in CEMGMsFEM, in general, do not have physical meanings since they are coordinates in the multiscale space. In order to obtain an upscaled model, we design multiscale basis functions such that the degrees of freedom represent the average of the solutions. As a result, the coarsegrid model is an equation for the solutions averaged over each continua. In this paper, we extend this idea to nonlinear equations.
To extend the concept of NLMC to nonlinear equations, we first identify macroscopic quantities for each coarsegrid block. These variables are typically found via local spectral decomposition and represent the features that can not be localized (similar to multicontina variables [6, 46, 53, 50]). Next, we consider local problems formulated in the oversampled regions with constraints. These local problems allow identifying the downscaling map from average macroscopic quantities to the finegrid variables. By imposing the constraints for each continuum variable via the source term, we define effective fluxes and the homogenized equation. Using the local solutions in the oversampled regions with constraints allows localizing the global downscaled map, which provides an accurate representation of the solution; however, it is expensive as it involves solving the global problem. By using the constraints in the oversampled regions, we can guarantee the proximity between the global and local downscaled maps for a given set of oversampled constraints.
The resulting homogenized equation significantly differs from standard homogenization. First, there are several variables per coarse block, which represent each continuum. Secondly, the local problems are formulated in oversampled regions with constraints. Finally, the nonlinear homogenized fluxes depend on all averages in oversampled regions, which bring nonlocal behavior for the equation. These ingredients are needed to perform upscaling in the absence of scale separation and high contrast.
In summary, our nonlinear nonlocal multicontinua approach has the following steps.

For each coarsegrid block, identify coarsegrid variables (continua) and associated quantities. This is typically done via some local spectral problems and describes the quantities that can not be localized.

Define local problems in oversampled regions with constraints. The constraints are given for each continua variables over all coarsegrid blocks. Local problems use specific source terms and boundary conditions, which allow localizing the global downscaled map and identifying effective fluxes. This step gives a downscaling from average continua to the finegrid solution.

The formulation of the coarsegrid problem. We seek the coarsegrid variables such that the downscaled finegrid solution approximately solves the global problem in a weak sense. The weak sense is defined via specific test functions, which are piecewise constants in each continua. This provides an upscaled model.
We consider two types of methods. First, we call linear interpolation and the second is nonlinear interpolation. In the first approach, we seek the solution in the form
where are multiscale basis functions for coarse cell and for continua defined in the oversampled region. This approach is simpler as one uses linear approximation of the solution. However, because of the nonlinearity, the nonlinear approximation is needed (see [29, 31, 30]). In [30], the authors propose such approach for pseudomonotone equations (similar to homogenization). In our problem, we present a nonlinear interpolation (we refer as nonlinear nonlocal multicontinua). In this approach, the solution is sought as a nonlinear map, which approximated from neighboring cells in a nonlinear fashion, which is defined via local problems. The local problem provides a nonlinear map from coarsegrid macroscopic variables to the finegrid solution. This mapping is used to construct macroscopic equations. We use nonlocality and multicontinua to address the cases without scale separation and high contrast.
As one of the numerical examples, we consider twophase flow and transport, though our approach can be applied to a wide range of problems. Twophase flow and transport is one of challenging problems, where many attempts are made to study it. Typical approaches include pseudorelative permeability approach, also known as multiphase upscaling, where the relative permeabilities are computed based on local twophase flow simulations. It is known that these approaches are process dependent. Our approach provides a novel upscaling, which shows that each coarsegrid block needs to contain several average pressures and saturations and, moreover, these coarsegrid relative permeabilities are nonlocal and depend on saturations and pressures of neighboring cells. This model provides a new way to look at twophase flow equations and can be further used in different modeling purposes.
The paper is organized as follows. In the next section, Section 2, we present some preliminaries and discuss main concepts used for linear problems. Section 3 is devoted to a general concept of nonlinear upscaling and this is discussed on several examples. In Section 4, we present the linear interpolation based approach and numerical results for multiphase flow and transport. Though the paper’s main focus is on nonlinear interpolation, the linear interpolation approach is practical for many applications. In Section 5, we present the nonlinear approach for multiphase flow and corresponding numerical results.
2 Preliminaries
We consider a general nonlinear system of the form
(1) 
In general, is a vectorvalued function and is the source term. is assumed to be heterogeneous in space (and time, in general), which are resolved on the fine grid (see Figure 1 for fine and coarse grid illustrations). Our objective is to derive coarsegrid equations. First, we discuss some examples for (1).
Example 1. Nonlinear (pseudomonotone) parabolic equations, where .
Example 2. Hyperbolic equations, where and is a scalar function and is a vector valued function. For example, and , as in Darcy’s flow.
Example 3. HamiltonJacobi equations, where and is a scalar function.
Example 4. One of our objective in this paper is to address the upscaling of multiphase flow, which is a challenging problem. In this case, the model equations have the following form. The twophase flow equations are derived by writing Darcy’s flow for each water () and oil () phases and the mass conservation as follows
(2) 
Here, is the Darcy velocity of the phase , , is the saturation of the phase , and is the relative mobility of the phase . By denoting the saturation of the water phase via , we can write the equations in the following way
(3) 
The special case of the twophase flow consists of singlephase flow
(4) 
2.1 Overview of NLMC for linear problems
In this section, we will give a brief overview of the NLMC method for linear problems [18]. Our goal is to summarize the key ideas in linear NLMC and motivate the ideas in the extension of it to the nonlinear NLMC in the next section. To be specific, we will consider a model parabolic equation with a heterogeneous coefficient, namely,
(5) 
with appropriate initial and boundary condition, where is the heterogeneous field, is a given source, is the physical domain and is a fixed time.
The NLMC upscaled system is defined on a coarse mesh, , of the domain . We write , where denotes the th coarse element and denotes the number of coarse elements in . For each coarse element , we will identify multiple continua corresponding to various solution features. This can be done via a local spectral problem or a suitable weight function. The upscaled parameter is defined using multiscale basis functions. For each coarse element and each continuum within , we will construct a multiscale basis function whose support is an oversampled region , which is obtained by enlarging the coarse block by a few coarse grid layers. See Figure 1 for an illustration of coarse grid and oversample region. In particular, a structured coarse grid is shown with boundaries of coarse elements are denoted blue. A coarse cell is denoted red and its oversampled region obtained by enlarging by two coarse grid layers is denoted green.
Now we will specify the definition of continuum. For each coarse block , we will identify a set of continua which are represented by a set of auxiliary basis functions , where denotes the th continuum. There are multiple ways to construct these functions . One way is to follow the idea proposed in CEMGMsFEM [19]. In this framework, the auxiliary basis functions are obtained as the dominant eigenfunctions of a local spectral problem defined on . These eigenfunctions can capture the heterogeneities and the contrast of the medium. The resulting solution of the upscaled system corresponds to the moments of the true solution with respect to these eigenfunctions. Another way is to follow the framework in the original NLMC method [18], designed for flows in fractured media, which can be easily modified for general heterogeneous media. In this approach, one identifies explicit information of fracture networks. The auxiliary basis functions are piecewise constant functions, namely, they equal one within one fracture network and zero otherwise. The resulting solution of the upscaled system corresponds to the average of the solution on fracture networks, which can be regarded as one of the continua. Finally, one can generalize the previous approach to construct other auxiliary basis functions. In particular, we can define to be characteristic functions of some regions . These regions can be chosen to reflect various properties of the medium. For instances, one can choose to be the region in which the medium has a certain range of values. We will consider this last choice in this paper.
Once the auxiliary basis functions are specified, we can construct the required basis functions. The idea generalizes the original energy minimization framework in CEMGMsFEM. Consider a given coarse element and a given continuum within . We will use the corresponding auxiliary basis function to construct our required multiscale basis function by solving a problem in an oversampled region . Specifically, we find and such that
(6)  
where denotes the standard delta function and is the space spanned by auxiliary basis functions. We remark the the function serves as a Lagrange multiplier for the constraints in the second equation of (6). We also remark that the basis function has mean value one on the th continuum within and has mean value zero in all other continua in all coarse elements within . In practice, the above system (6) is solved in using a fine mesh, which is typically a refinement of the coarse grid. See Figure 1 for an illustration.
Now, we can derive the NLMC upscaled system. As an illustration of the concept, we consider a forward Euler method for the time discretization of (5). At the th time step, the solution vector is denoted by , where each component of represents the average of the solution on a continuum within a coarse element. We note that the size of this vector is , where is the number of continua in . The upscaled stiffness matrix is defined as
(7) 
and the upscaled mass matrix is defined as
(8) 
Finally, the NLMC system is written as
(9) 
where the components of the vector are defined as and is the time at the th time step. We remark that the nonlocal connections of the continua are coupled by the matrices and . We also remark that the local computation in (6) results from a spatial decay property of the multiscale basis function, see [19, 20, 16] for the theoretical foundation.
3 Nonlinear nonlocal multicontinua model
3.1 General concept
We will first present some general concept of our nonlinear NLMC. We consider the following nonlinear problem
(10) 
where has a multiscale dependence with respect to space (and time, in general).
To formulate the coarsegrid equations in the time interval , we first introduce coarsegrid variables , where is the coarsegrid block, is a continuum representing the coarsegrid variables, and is the time step (cf. [18]). As we noted that the continuum plays a role of a macroscopic variable, which can not be localized. For each coarsegrid block , we need several coarsegrid variables, which will be indexed by . In this paper, we will consider two types of interpolation for multiscale degrees of freedom. The first approach is called linear interpolation, which constructs approximate solution at the time as
where are multiscale basis functions, which are defined via local constrained problems formulated in the oversampled regions. These basis functions are possibly supported in the oversampled regions and the resulting coarsegrid equations will be nonlinear and nonlocal when projected to the appropriate test spaces. The values of are computed via the variational formulation using an explicit or implicit discretization of time
(11) 
where or depending whether explicit or implicit discretization is used, and denotes test functions. Here, denotes usual inner product. As for the test space, one can use the multiscale space spanned by or the auxiliary basis functions . This linear approach is simpler to use; however, it ignores the nonlinear interpolation, which is important for nonlinear homogenization and numerical homogenization [51, 30].
Next, we describe the nonlinear approach, which we refer as nonlinear nonlocal multicontinuum approach. In this case, the solution is sought as a nonlinear interpolation of the degrees of freedom defined in the oversampled region. More precisely, we compute a downscale function using the given solution values by solving a constrained problem in the oversampled regions. Mathematically, this downscale function can be written as a function of all values , namely,
where is a nonlinear map and is a vector containing all values . The nonlinear map is computed by solving the local problem (cf. (14)), and the function is glued together from local downscaled maps. Once this map is defined, we seek all coarsegrid macroscopic variables such that the downscaled finegrid solution solves the global problem in a variational setting. The test functions are defined for each degrees of freedom in the form of piecewise constants or piecewise functions, in for example a finite volume or Petrov Galerkin setting. More precisely, the coarsegrid system has the form
(12) 
for all suitable test functions . Here, or can be used for explicit or implicit discretization. The test functions are chosen to be piecewise polynomials or auxiliary basis functions . We remark that the dimension of the test space is chosen to be the dimension of the coarsegrid macroscopic variables.
3.2 Nonlinear nonlocal multicontinuum approach
In this section, we give some details of our nonlinear nonlocal multicontinuum approach, in general, and then present some examples. The illustration is given in Figure 2. In later sections, we present a detailed application to twophase flow equations.
Step 1. Defining coarsegrid variables. The first step involves defining macroscopic variables. In our cases, we will define them by , where is the coarsegrid block, is the continua, and is the time step. These variables in our applications will be defined by prescribing averages to subregions, which can have a complex shapes (for example, fracture regions). The equation for , in general, will have a form
(13) 
where or for explicit or implicit discretizations. The operator is determined respect to the given continuum and coarse element , and contains information from the source term and the time step size . In (13), we use to denote the vector consisting of all macroscopic quantities for all and . The operator is obtained by solving local problems on an oversampled region corresponding to . To define (13), we next discuss local problems.
Step 2. Local solves for generic constraints. The computation of requires solutions of constrained local nonlinear problems. Here we present a general formulation. We consider to be the index for a coarse element or a coarse neighborhood , which is defined for a coarse node by
The choice of or depends on the global coarse grid discretization. For example, with finite volume or Petrov Galerkin formulation, we choose , while for continuous Galerkin formulation, we choose . In this part, we present the steps using . The required local nonlinear problem will be solved on , which is an oversampled region obtained by enlarging a few coarse grid layers. We let be a set of scalar values, where denotes the th coarse element in the oversampled region and denotes the th continuum within . The local problem is solved by finding a function given by
(14) 
with constraints
Here, is an indicator function for the region defined for the continua within coarse block . We remark that the values play the role of Lagrange multipliers for the constraints. We notice that depends on all constraints in the oversampled region. In general, one can also impose constraints on the gradients of and the constraint function can have a complex form. We remark that the problem (14) requires a boundary condition, which is problem dependent. We will discuss it later. Our main message is that the local problems with constraints are solved in oversampled region and involves several coarsegrid variables for each coarse block. For the precise formulation, the values will be chosen as the macroscopic variables.
Discussion on localization.
Next, we discuss the idea behind the localization principle stated above (see Figure 3 for illustration). For this reason, we first define the global downscaling operator, , which solves the global problem with constraints
(15) 
with constraints
The boundary conditions are taken as the original global boundary conditions. For our localization, we desire that the local downscaled solution, in (target coarse block) defined in (14), is approximately the same as restricted to for the same constraint values as the global problem in . I.e.,
where values coincide in . This property allows solving local nonlinear problems instead of the global problem without sacrificing the accuracy and defining correct upscaled coefficients. This desired property can be shown for monotone elliptic equations using zero Dirichlet boundary conditions on . For more complex nonlinear system, one needs to also use appropriate boundary conditions on to achieve this property.
Step 3. Defining coarsegrid model. The upscaled flux is defined by substituting the downscaled global solution corresponding to the constraints given by and multiplying it by test functions. For this procedure, one first defines a global finegrid downscaled field that is constrained. The equation (14) defines local finegrid fields constrained to macroscopic variables in oversampled regions. In general, one needs to “glue” together a global downscaled solution, which approximates the global problem in a weak sense. A simple approach to glue together is to use partition of unity functions, for each region . Then,
In some discretization scenarios (e.g., finite volume, Discontinuous Galerkin,…), it is not necessary to glue together in order to obtain a global finegrid field. In the proposed NLMC approach, our coarsegrid fomrulation uses “finite volume” type approximation and we will not need to glue local approximations. In this case,
(16) 
where is a coarse edge for the coarse grid. In general, we can write it as in all local finegrid blocks, and the following variational formulation
(17) 
where are test functions. In this setup,
The above equation can formally be thought as
The time discretization of Equation (17) can be written as
(18) 
where or for explicit or implicit discretizations. In general, the downscaling operator has the following property
(19) 
which simplifies the computations. In the latter situation, the mass matrix is diagonal. In addition,
(20) 
where is a finescale field and is the corresponding coarsegrid average. The equation (20) states that if we use the average of the finescale field to reconstruct it, the resulting approximation remains close to the original finescale field. We remark that the operator in (13) is obtained using (18). More precisely, we have
by using a test function corresponding to region and continuum .
Next, we present some more concrete constructions for some cases, listed in Section 2.
Example 1. In this example, we consider a class of nonlinear (pseudomonotone) parabolic equations with , where is a heterogeneous function in , and in general, also in and . We will consider a continuous Galerkin discretization in space on a coarse grid and explicit Euler discretization in time. Let be a set of macroscopic variables at the time . In particular, denotes the mean value of the solution on continuum within coarse element at time time . On an oversampled region , we find a function on the fine grid such that
subject to the following constraints
We remark that plays the role of Lagrange multiplier. The above problem is solved using the Dirichlet boundary condition on . We can set equals to zero on , and this choice of motivated by the decay property of multiscale basis functions in CEMGMsFEM. Then we can define a global fine scale function by
where are suitable partition of unity functions, where we notice that the values of within are not used. Using the global downscale function and suitable test functions, we can derive a coarse grid scheme. We note that the upscaled model has the form (17).
Example 2. In this example, we consider a class of hyperbolic equations, where and is a scalar function and is a vector valued function. For example, and , as in Darcy’s flow. We will consider a finite volume type discretization in space on a coarse grid and explicit Euler discretization in time. Let be a set of macroscopic variables at the time . Same as above, denotes the mean value of the solution on continuum within coarse element at time time . On an oversampled region , we find a function on the fine grid such that
subject to the following constraints
We remark that plays the role of Lagrange multiplier. The above problem is solved using the standard inflow boundary condition on the inflow part of . We can set equals to zero to the value chosen by upwinding. Then we can define a global fine scale function by
where are suitable partition of unity functions for the overlapping partition . We remark that this global downscale function preserves the mean values, namely,
Using the global downscale function and suitable test functions, we can derive a coarse grid scheme. For the test functions, we use , for all . Hence, we obtain
(21) 
We note that this upscaled model has the form (17). We also note that the number of unknowns is the same of number of equations.
Example 3. For this case, the local problem can be formulated as follows. On an oversampled region , we find a function on the fine grid such that
subject to the following constraints
Similar to Example 2, the upscaled model has the form (17) and (21).
In summary, the local problems involve original local problems with constraints formulated in the oversampled regions. The source term constants represent the homogenized fluxes and their dependence on averages of solutions and gradients are the functional form of the equations.
3.2.1 RVEbased extension of the method
The proposed concept can also be used for problems with scale separation. A typical example is a fractured media, where the fracture distributions are periodic or posess some scale separation. In this case, we do not use the oversampling based local problems and simply use local problems in RVE. To be more precise, Step 1 remains the same as before, which involves identifying local multicontinua variables in each coarsegrid block.
In Step 2, we solve the local problem (14) in RVE subject to the constraints. The main difference is that we use only the constrained in the target coarsegrid block. Thus, there are fewer constraints and the problem is localized to the target coarse block.
4 Linear Approach
In this section, we will focus on the linear approach (c.f. Section 3), and present some numerical results.
4.1 Linear transport
In this section, we will present some numerical results for the linear transport equation with a given velocity. More precisely, we consider
where is a given divergencefree velocity field with on . We take , and the velocity field is shown in Figure 4. The coarse grid size , and we use a single continuum model. In Table 1, we present a error comparison between our upscaling method and the standard finite volume method.
The derivation of the upscaled system follows the ideas in Section 2.1. For each coarse element , we solve the following system in an oversampled region :
(22) 
to obtain a basis function . The above system is equipped with the constraints
(23) 
We remark that is a piecewise constant function, and plays the role of Lagrange multiplier. We notice that the upscaled system has the form (9). In practice, one can add artificial diffusion in (22) in order to obtain a stable numerical solution, and use zero Dirichlet boundary condition. We remark that one can use other boundary conditions, see Section 3.2.
By using the standard finite volume method, the relative error is at the final time . From Table 1, we see that our upscaling method provides much better approximations, where ” layer” standards for the number of layers used in the oversampling domains. When layer equals , it means that the oversampling domain is the whole physical domain. In Figure 5, we show the snapshots of the solution at . From this figure, we observe that the solution computed with CEM provides a good approximation, particularly, in saturated regions.
#layer  Errors 

5  11.57% 
7  7.73% 
9  7.22% 
6.91% 
4.2 Single phase flow
In this section, we apply our method to a single phase flow problem which is described as follows:
(24)  
(25)  
(26) 
where we assume is a function independent of time such that we can compute the velocity in the initial time step and use the velocity to construct the basis function for saturation .
For solving single phase flow problem, we can separate the computation process into two parts. Since Equations (25) and (26) are independent of saturation and the given source term is independent of time, we know that the velocity is independent of time. Therefore, we can decouple system of equation into two parts which are the Equations (25)(26) and the Equation (24). The first part of the computation is computing the velocity and pressure by solving Equations (25) and (26) by Mixed Generalized Multiscale Finite Element Methods(MGMsFEM) [11]. The second part of the computation is computing the numerical solution for the saturation by solving Equation (24). Given a numerical velocity , Equation (24) is a standard linear transport equation. We can use the method discussed in Section 4.1 to construct the multiscale basis function and then compute the numerical solution for saturation. Next, we will discuss the computation of the saturation.
We will consider a multicontinuum version of the method. Assume that an approximation of the velocity has been computed by the MGMsFEM. We derive the upscaled system following the ideas in Section 2.1. For each coarse element and a continuum within , we solve the following system in an oversampled region :
(27) 
to obtain a basis function . The above system is equipped with the constraints
(28) 
We remark that is a piecewise constant function, and plays the role of Lagrange multiplier. We notice that the upscaled system has the form (9).
Now, we present some numerical examples. We take , and the permeability field is shown in Figure 6. The coarse grid size . We will use a dual continuum model to solve this problem. For each coarse element , we define two continua by and . In Table 2, we present a error comparison at times and between our method and the finite volume method for the saturation equation (26) with piecewise constant test functions in each continuum. From the table, we observe the performance of our scheme for various choices of oversampling layers. For all these cases, we observe that our scheme performs better than the usual (dual continuum) finite volume scheme whose error is and at the observation times and respectively. We notice that the relative error in the velocity is . In Figure 7  8, we present the comparisons between the the averaged finegrid solution and the solution obtained our proposed method. From these figures, we observe that our proposed method can capture the solution features accurately.
#layer \ time  

4  12.83%  14.57% 
6  7.33%  6.93% 
8  5.89%  5.51% 
5.69%  5.24% 
4.3 Two phase flow
In this section, we apply our upscaling method to two phase flow problem which is described as follows:
(29)  
(30)  
(31)  
(32) 
In the process of basis construction, we will use a similar approach as solving single phase flow problem. We notice that the exact velocity depends on the saturation and hence depends on time. To reduce the computational burden, we use the initial total mobility as the reference mobility to construct the multiscale basis functions. That is, we will construct the multiscale basis functions for approximating the exact velocity field based on the following equations,
(33) 
(34) 
Next, using the total velocity , we will construct the basis functions for approximating the exact saturation based on the following transport equation,
(35) 
Hence, the construction for the basis functions is separated into two parts. The first part is computing an approximating velocity and pressure by solving Equations (33) and (34) by Mixed Generalized Multiscale Finite Element Methods (MGMsFEM). The multiscale basis functions constructed during the processing are used to span the finite element space , which is used to approximate the exact velocity. The second part is constructing the saturation basis functions by the method discussed in Section 3 with a fixed velocity . This part is also similar to the single phase flow.
To compute the coarse grid solution, we will use IMplicit Pressure Explicit Saturation (IMPES) scheme. Given the saturation and the total velocity at the previous time step, we will compute the total velocity at the next time step by solving
where is the space for pressure. Then, we will use the velocity to compute the saturation at the next time step as before (similar to (27) and (28)).
We next present some numerical examples. We take . The permeability field is shown in Figure 9, and we use
with , and . The coarse grid size . We are using the dual continuum model to solve this problem. For each coarse element , we define two continua by and . In Table 3, we present a error comparison between our method and the finite volume method for the saturation equation with piecewise constant test functions in each continuum. For this experiment, we observed that the upscaling method works in general better than the finite volume method as before.
#layer \ time  

4  15.82%  14.94% 
6  7.03%  6.63% 
8  4.11%  4.33% 
3.52%  4.11% 
Next we consider another test case with the medium shown in Figure 10, which is a SPE benchmark test case. The coarse grid size is chosen as . We are again using the dual continuum model to solve this problem. For each coarse element , we define and . In Figure 11, we present the snapshots of the averaged fine solution, finite volume solution and our upscaled solution at the observation time . We observe that our method is able to capture the dynamics of the solution, while the finite volume method does not produce a good solution. In Table 4, we present the errors of our approximation and observe again good accuracy.
#layer \ time  

6  7.19%  11.12% 
6.70%  10.16% 
5 Nonlinear Approach
In this section, we will present a numerical test to show the performance of the nonlinear approach (c.f. Section 3). To to do, we consider the following single phase flow problem
where
We assume is a function independent of time such that we can compute the velocity in the initial time step and use the velocity to construct the basis function for saturation.
The derivation of the upscaled system follows the ideas in Section 3.2. Assume that an approximation of the velocity has been computed. Let be a set of upscaled values for the continuum in the coarse region at the time . For each coarse element , we solve the following system in an oversampled region :
(36) 
to obtain a local downscale function . The above system is equipped with the constraints
(37) 
We remark that is a piecewise constant function, and plays the role of Lagrange multiplier. Then we define a global downscale field by
where is a set of partition of unity functions corresponding to the sets . Then we apply a finite volume scheme to the saturation equation as follows
(38) 
We next present some numerical results. We take , and the permeability field is shown in Figure 12. The coarse grid size . We are using a dual continuum model to solve this problem. For each , we define and . In Table 5, we present a error comparison between our method and the finite volume method for the saturation equation with piecewise constant test functions in each continua. We consider three choices of , and present the relative errors at times and . From the results, we observe that our method is able to compute more accurate numerical solutions.
\ time  

2  9.91%  11.85% 
3  7.54%  12.60% 
5  10.07%  14.89% 
\ time  

2  28.70%  23.61% 
3  31.32%  26.27% 
5  32.98%  29.19% 
We remark that our proposed upscaling method can be applied to other nonlinear systems, which will be investigated in more detail in a forthcoming paper.
6 Conclusions
In this paper, we study and develop a general framework for coarsegrid modeling of nonlinear PDEs. The finegrid problem is described by nonlinear PDEs including twophase flow and transport model. The main concept of our upscaled model is three fold, (1) We determine macroscopic quantities in each coarsegrid block, which can vary among the coarsegrid blocks. (2) We define appropriate local problems with constraints in oversampled regions which are solvable. The constraints are given by macroscopic quantities. These local problems allow downscaling from macroscopic quantities to each coarsegrid block. (3) We formulate the coarsegrid problem such that the downscaled solution over the entire domain solves the global problem in a weak sense.
The formulation of local problems allow easily defining macroscopic fluxes. The local problems are defined via forcing terms. We also discuss the use of multiscale basis function in a linear fashion for solving nonlinear PDEs. Our main example includes twophase flow and transport and its simplifications. The model problem is described by a coupled PDEs. We consider various permeability fields and show that one can achieve an accurate approximation of the solution.
Acknowledgements
EC’s work is partially supported by Hong Kong RGC General Research Fund (Project 14304217) and CUHK Direct Grant for Research 201718. EC would like to thank the support of the ESI for attending the program Numerical Analysis of Complex PDE Models in the Sciences.
References
 [1] Assyr Abdulle and Yun Bai. Adaptive reduced basis finite element heterogeneous multiscale method. Comput. Methods Appl. Mech. Engrg., 257:203–220, 2013.
 [2] G. Allaire and R. Brizzi. A multiscale finite element method for numerical homogenization. SIAM J. Multiscale Modeling and Simulation, 4(3):790–812, 2005.
 [3] T. Arbogast. Implementation of a locally conservative numerical subgrid upscaling scheme for twophase Darcy flow. Comput. Geosci, 6:453–481, 2002.
 [4] T. Arbogast, G. Pencheva, M.F. Wheeler, and I. Yotov. A multiscale mortar mixed finite element method. SIAM J. Multiscale Modeling and Simulation, 6(1):319–346, 2007.
 [5] T. Arbogast, G. Pencheva, M.F. Wheeler, and I. Yotov. A multiscale mortar mixed finite element method. Multiscale Model. Simul., 6(1):319–346, 2007.
 [6] GI Barenblatt, Iu P Zheltov, and IN Kochina. Basic concepts in the theory of seepage of homogeneous liquids in fissured rocks [strata]. Journal of applied mathematics and mechanics, 24(5):1286–1303, 1960.
 [7] J.W. Barker and S. Thibeau. A critical review of the use of pseudorelative permeabilities for upscaling. SPE Reservoir Eng., 12:138–143, 1997.
 [8] Donald L Brown and Daniel Peterseim. A multiscale method for porous microstructures. arXiv preprint arXiv:1411.1944, 2014.
 [9] Y. Chen, L. Durlofsky, M. Gerritsen, and X. Wen. A coupled localglobal upscaling approach for simulating flow in highly heterogeneous formations. Advances in Water Resources, 26:1041–1060, 2003.
 [10] E. Chung, Y. Efendiev, and S. Fu. Generalized multiscale finite element method for elasticity equations. International Journal on Geomathematics, 5(2):225–254, 2014.
 [11] E. Chung, Y. Efendiev, and C. Lee. Mixed generalized multiscale finite element methods and applications. SIAM Multicale Model. Simul., 13:338–366, 2014.
 [12] E. Chung, Y. Efendiev, and W. T. Leung. Generalized multiscale finite element method for wave propagation in heterogeneous media. SIAM Multicale Model. Simul., 12:1691–1721, 2014.
 [13] E. Chung and W. T. Leung. A subgrid structure enhanced discontinuous galerkin method for multiscale diffusion and convectiondiffusion problems. Communications in Computational Physics, 14:370–392, 2013.
 [14] E. T. Chung, Y. Efendiev, W.T. Leung, M. Vasilyeva, and Y. Wang. Online adaptive local multiscale model reduction for heterogeneous problems in perforated domains. Applicable Analysis, 96(12):2002–2031, 2017.
 [15] E. T. Chung, Y. Efendiev, and G. Li. An adaptive GMsFEM for high contrast flow problems. J. Comput. Phys., 273:54–76, 2014.
 [16] Eric Chung, Yalchin Efendiev, and Wing Tat Leung. Constraint energy minimizing generalized multiscale finite element method in the mixed formulation. Computational Geosciences, 22(3):677–693, 2018.
 [17] Eric Chung, Maria Vasilyeva, and Yating Wang. A conservative local multiscale model reduction technique for stokes flows in heterogeneous perforated domains. Journal of Computational and Applied Mathematics, 321:389–405, 2017.
 [18] Eric T Chung, Efendiev, Wing Tat Leung, Maria Vasilyeva, and Yating Wang. Nonlocal multicontinua upscaling for flows in heterogeneous fractured media. arXiv preprint arXiv:1708.08379, 2018.
 [19] Eric T Chung, Yalchin Efendiev, and Wing Tat Leung. Constraint energy minimizing generalized multiscale finite element method. Computer Methods in Applied Mechanics and Engineering, 339:298–319, 2018.
 [20] Eric T Chung, Yalchin Efendiev, and Wing Tat Leung. Fast online generalized multiscale finite element method using constraint energy minimization. Journal of Computational Physics, 355:450–463, 2018.
 [21] Martin Drohmann, Bernard Haasdonk, and Mario Ohlberger. Reduced basis approximation for nonlinear parametrized evolution equations based on empirical operator interpolation. SIAM J. Sci. Comput., 34(2):A937–A969, 2012.
 [22] L.J. Durlofsky. Numerical calculation of equivalent grid block permeability tensors for heterogeneous porous media. Water Resour. Res., 27:699–708, 1991.
 [23] W. E and B. Engquist. Heterogeneous multiscale methods. Comm. Math. Sci., 1(1):87–132, 2003.
 [24] Y. Efendiev and L.J. Durlofsky. Numerical modeling of subgrid heterogeneity in two phase flow simulations. Water Resour. Res., 38(8):1128, 2002.
 [25] Y. Efendiev and L.J. Durlofsky. A generalized convectiondiffusion model for subgrid transport in porous media. SIAM J. Multiscale Modeling and Simulation, 1(3):504–526, 2003.
 [26] Y. Efendiev, J. Galvis, and T. Y. Hou. Generalized multiscale finite element methods (gmsfem). Journal of Computational Physics, 251:116–135, 2013.
 [27] Y. Efendiev, J. Galvis, and X.H. Wu. Multiscale finite element methods for highcontrast problems using local spectral basis functions. Journal of Computational Physics, 230:937–955, 2011.
 [28] Y. Efendiev and T. Hou. Multiscale Finite Element Methods: Theory and Applications, volume 4 of Surveys and Tutorials in the Applied Mathematical Sciences. Springer, New York, 2009.
 [29] Y. Efendiev and A. Pankov. Numerical homogenization of monotone elliptic operators. SIAM J. Multiscale Modeling and Simulation, 2(1):62–79, 2003.
 [30] Y. Efendiev and A. Pankov. Numerical homogenization of nonlinear random parabolic operators. SIAM J. Multiscale Modeling and Simulation, 2(2):237–268, 2004.
 [31] Y. Efendiev and A. Pankov. Homogenization of nonlinear random parabolic operators. Advances in Differential Equations, 10(11):1235–1260, 2005.
 [32] DA Fafalis, SP Filopoulos, and GJ Tsamasphyros. On the capability of generalized continuum theories to capture dispersion characteristics at the atomic scale. European Journal of MechanicsA/Solids, 36:25–37, 2012.
 [33] Dimitrios Fafalis and Jacob Fish. Computational continua for linear elastic heterogeneous solids on unstructured finite element meshes. International Journal for Numerical Methods in Engineering, 115(4):501–530, 2018.
 [34] Jacob Fish. Practical multiscaling. John Wiley & Sons, 2013.
 [35] Jacob Fish and Wen Chen. Space–time multiscale model for wave propagation in heterogeneous media. Computer Methods in applied mechanics and engineering, 193(45):4837–4856, 2004.
 [36] Jacob Fish and Rong Fan. Mathematical homogenization of nonperiodic heterogeneous media subjected to large deformation transient loading. International Journal for numerical methods in engineering, 76(7):1044–1064, 2008.
 [37] Jacob Fish, Vasilina Filonova, and Dimitrios Fafalis. Computational continua revisited. International Journal for Numerical Methods in Engineering, 102(34):332–378, 2015.
 [38] Jacob Fish, Vasilina Filonova, and Zheng Yuan. Reduced order computational continua. Computer Methods in Applied Mechanics and Engineering, 221:104–116, 2012.
 [39] Jacob Fish and Sergey Kuznetsov. Computational continua. International Journal for Numerical Methods in Engineering, 84(7):774–802, 2010.
 [40] Jacob Fish, Kamlun Shek, Muralidharan Pandheeradi, and Mark S Shephard. Computational plasticity for composite structures based on mathematical homogenization: Theory and practice. Computer Methods in Applied Mechanics and Engineering, 148(12):53–73, 1997.
 [41] Jacob Fish and Zheng Yuan. Multiscale enrichment based on partition of unity. International Journal for Numerical Methods in Engineering, 62(10):1341–1359, 2005.
 [42] Jacob Fish and Zheng Yuan. Multiscale enrichment based on partition of unity for nonperiodic fields and nonlinear problems. Computational Mechanics, 40(2):249–259, 2007.
 [43] Patrick Henning and Mario Ohlberger. The heterogeneous multiscale finite element method for elliptic homogenization problems in perforated domains. Numerische Mathematik, 113(4):601–629, 2009.
 [44] L. Holden and B.F. Nielsen. Global upscaling of permeability in heterogeneous reservoirs: the Output Least Squares (OLS method. Transport in Porous Media, 40:115–143, 2000.
 [45] J.R. Kyte and D.W. Berry. New pseudofunctions to control numerical dispersion. Society of Petroleum Engineers Journal, 15(4):269–276, 1975.
 [46] Seong H Lee, MF Lough, and CL Jensen. Hierarchical modeling of flow in naturally fractured formations with multiple length scales. Water resources research, 37(3):443–455, 2001.
 [47] AnaMaria Matache and Christoph Schwab. Twoscale fem for homogenization problems. ESAIM: Mathematical Modelling and Numerical Analysis, 36(04):537–572, 2002.
 [48] Caglar Oskay and Jacob Fish. Eigendeformationbased reduced order homogenization for failure analysis of heterogeneous materials. Computer Methods in Applied Mechanics and Engineering, 196(7):1216–1243, 2007.
 [49] H. Owhadi and L. Zhang. Metricbased upscaling. Comm. Pure. Appl. Math., 60:675–723, 2007.
 [50] GP Panasenko. Multicontinuum wave propagation in a laminated beam with contrasting stiffness and density of layers. Journal of Mathematical Sciences, pages 1–13, 2018.
 [51] A. Pankov. convergence and homogenization of nonlinear partial differential operators. Kluwer Academic Publishers, Dordrecht, 1997.
 [52] M. Peszyńska, M. Wheeler, and I. Yotov. Mortar upscaling for multiphase flow in porous media. Comput. Geosci., 6(1):73–100, 2002.
 [53] JE Warren, P Jj Root, et al. The behavior of naturally fractured reservoirs. 1963.
 [54] X.H. Wu, Y. Efendiev, and T.Y. Hou. Analysis of upscaling absolute permeability. Discrete and Continuous Dynamical Systems, Series B., 2:158–204, 2002.
 [55] Zheng Yuan and Jacob Fish. Multiple scale eigendeformationbased reduced order homogenization. Computer Methods in Applied Mechanics and Engineering, 198(2126):2016–2038, 2009.