Space-Time Nonlinear Upscaling Framework Using Non-local Multi-continuum Approach

TL;DR

This paper introduces a space-time nonlinear upscaling framework for porous media that handles multiscale coefficients without scale separation, utilizing nonlocal multi-continuum concepts and machine learning for efficient solutions.

Contribution

The paper extends previous results by developing a nonlinear nonlocal multi-continuum upscaling method applicable to complex porous media problems without scale separation.

Findings

Effective upscaling for nonlinear porous media problems.

Use of machine learning to identify complex solution maps.

Numerical validation on two-phase flow and transport applications.

Abstract

In this paper, we develop a space-time upscaling framework that can be used for many challenging porous media applications without scale separation and high contrast. Our main focus is on nonlinear differential equations with multiscale coefficients. The framework is built on nonlinear nonlocal multi-continuum upscaling concept and significantly extends the results in the proceeding paper. Our approach starts with a coarse space-time partition and identifies test functions for each partition, which plays a role of multi-continua. The test functions are defined via optimization and play a crucial role in nonlinear upscaling. In the second stage, we solve nonlinear local problems in oversampled regions with some constraints defined via test functions. These local solutions define a nonlinear map from macroscopic variables determined with the help of test functions to the fine-grid fields. This map can be thought as a downscaled map from macroscopic variables to the fine-grid solution. In the final stage, we seek macroscopic variables in the entire domain such that the downscaled field solves the global problem in a weak sense defined using the test functions. We present an analysis of our approach for an example nonlinear problem. Our unified framework plays an important role in designing various upscaled methods. Because local problems are directly related to the fine-grid problems, it simplifies the process of finding local solutions with appropriate constraints. Using machine learning (ML), we identify the complex map from macroscopic variables to fine-grid solution. We present numerical results for several porous media applications, including two-phase flow and transport.

Tables

Table 1. Table 1 : Learning performance of machine learning algorithm for Test 1 and Test 2

	MSE	RMSE (%)	MAE (%)
Test 1
$N{N}_1$	0.113	3.368	2.798
$N{N}_2$	0.029	1.725	1.587
$N{N}_3$	0.283	5.322	4.381
$N{N}_4$	0.048	2.196	2.443

Equations

$$-\nabla \cdot (\kappa \nabla u)=f,\text{ in }\Omega .$$

$${\int }_{K_i^{+}}\kappa \nabla {\psi }_i^j\cdot \nabla v+{\int }_{K_i^{+}}\tilde{\kappa }\mu v=0,\forall v\in H_0^1(K_i^{+}),$$

$${\int }_{K_{\ell }}\tilde{\kappa }{\psi }_i^j{\varphi }_m^{\ell }={\delta }_{j\ell }{\delta }_{im},\forall K_{\ell }\subset K_i^{+},$$

$$u_{ms}=\sum _{i=1}^N\sum _j{U}_i^j{\psi }_i^j.$$

$$a(u_{ms},\psi )=(g,\psi ),\forall \psi \in V_{ms}.$$

$$A_{T}U=F$$

$$(A_T{)}_{jm}^{(i,\ell )}=a({\psi }_j^i,{\psi }_m^{\ell }):={\int }_{\Omega }\kappa \nabla {\psi }_j^i\cdot \nabla {\psi }_m^{\ell },$$

$$(F{)}_i^{(j)}=(g,{\psi }_i^j).$$

$$M{U}_t+\nabla \cdot G(x,t,U)=g,$$

$$U_i^{(j)}=⟨⟨U,{\psi }_i^{(j)}⟩⟩$$

$$M{\varphi }_t+\nabla \cdot G(x,t,\varphi )=\mu ,\text{in }{K}_i^{+}.$$

$$I_{\varphi }({\psi }_i^{(j)}(x,t)).$$

$${\mathcal{F}}_i^{ms}:I_{\varphi }({\psi }_i^{(j)}(x,t))\to \varphi .$$

$$⟨⟨M{U}_t^{ms}+\nabla \cdot G(x,t,U^{ms}),{\psi }_i^{(j)}⟩⟩=⟨⟨g,{\psi }_i^{(j)}⟩⟩.$$

$$\begin{array}{r}M({\varphi }_{i,m}^{(j,l)}{)}_t+\nabla \cdot G(x,t,{\varphi }_{i,m}^{(j,l)})={\mu }_{i,m}^{(j,l)}\\ ⟨⟨{\varphi }_{i,m}^{(j,l)},{\psi }_s^{(r)}⟩⟩={\delta }_{lr}{\delta }_{ms}\end{array}$$

$${\psi }_i^{(j)}(x,t)=\sum _{m,l}{U}_m^{(l)}{\varphi }_{i,m}^{(j,l)}.$$

$$G(x,t,\sum _{m,l}{U}_m^{(l)}{\varphi }_{i,m}^{(j,l)})=\sum _{m,l}{U}_m^{(l)}G(x,t,{\varphi }_{i,m}^{(j,l)}).$$

$$⟨⟨M{\varphi }_t+\nabla \cdot G(x,t,\varphi ),{\psi }_i^{(j)}⟩⟩=⟨⟨g,{\psi }_i^{(j)}⟩⟩,\forall {\psi }_i^{(j)}.$$

$$\begin{array}{r}M({\mathcal{F}}^{ms}(U){)}_t+\nabla \cdot G(x,t,{\mathcal{F}}^{ms}(U))=\mu \\ ⟨⟨{\mathcal{F}}^{ms}(U),{\psi }_i^{(j)}⟩⟩=U_i^{(j)}\end{array}$$

$$\sum _{l,m}⟨⟨M({\mathcal{F}}^{ms}(U){)}_t+\nabla \cdot G(x,t,{\mathcal{F}}^{ms}(U)),{\psi }_i^{(j)}⟩⟩=⟨⟨g,{\psi }_i^{(j)}(x,t)⟩⟩\forall {\psi }_i^{(j)}.$$

$${\partial }_{t}u+L(u)=f,\text{in }\Omega \times (0,T]$$

$$U_j^{(n,i)}=s(u,{\psi }_j^{(n,i)})={\int }_0^T{\int }_{\Omega }\tilde{\kappa }{\psi }_j^{(n,j)}u$$

$${\int }_0^T{\int }_{\Omega }({\partial }_t+L)F_1(U)v-s(F_2(U),v)$$

$$s(F_1(U),{\psi }_j^{(n,i)})$$

$$s(F_2(U),{\psi }_j^{(n,i)})={\int }_0^T{\int }_{\Omega }f{\psi }_j^{(n,i)},\forall {\psi }_j^{(n,i)}\in V_{aux}.$$

$${\int }_{t_n^{-}}^{t_{n+1}}{\int }_{K_{+}^i}({\partial }_t+L)F_{loc,1}^{(n,i)}(U)v-s(F_{loc,2}^{(n,i)}(U),v)$$

$$s(F_{loc,1}^{(n,i)}(U),{\psi }_j^{(n,i)})$$

$$s(F_{ms,2}(U),{\psi }_j^{(n,i)})={\int }_0^T{\int }_{\Omega }f{\psi }_j^{(n,i)},\forall {\psi }_j^{(n,i)}\in V_{aux}$$

$$\nabla \cdot (\kappa (x,\nabla u))$$

$$u$$

$$A_{\Omega }(u,w)$$

$$A_{\omega }(u,w)={\int }_{\omega }\kappa (x,\nabla u)\cdot \nabla w.$$

Full text

Our approach starts with a coarse space-time partition and identifies test functions for each partition, which play a role of multi-continua. The test functions are defined via optimization and play a crucial role in nonlinear upscaling. In the second stage, we solve nonlinear local problems in oversampled regions with some constraints defined via test functions. These local solutions define a nonlinear map from macroscopic variables determined with the help of test functions to the fine-grid fields. This map can be thought as a downscaled map from macroscopic variables to the fine-grid solution. In the final stage, we seek macroscopic variables in the entire domain such that the downscaled field solves the global problem in a weak sense defined using the test functions. We present an analysis of our approach for an example nonlinear problem.

Our unified framework plays an important role in designing various upscaled methods. Because local problems are directly related to the fine-grid problems, it simplifies the process of finding local solutions with appropriate constraints [16]. Using machine learning (ML), we identify the complex map from macroscopic variables to fine-grid solution. We present numerical results for several porous media applications, including two-phase flow and transport.

1 Introduction

Many porous media models are nonlinear and deriving these nonlinear macroscopic equations rely on some assumptions. For example, the well-known two-phase flow and transport model assumes that the relative permeabilities are functions of local saturations [6]. Similarly, for unsaturated flows, the nonlinear relations between pressures and capillary curves use local relations. All these problems have space-time heterogeneities. Some rigorous upscaling tools are needed to generalize these models and understand the errors associated in these macroscopic models. This is one of our goals in this paper.

Many approaches are suggested for nonlinear upscaling in the past, e.g., [2, 26, 3, 25, 13, 7, 9, 22, 1, 20, 31, 32, 42, 40, 4, 37, 12, 15, 10, 44, 5, 11, 36, 41, 46]. For multi-phase flows, these techniques include permeability or transmissibility upscaling [21, 45, 8, 38] for single-phase flow and pseudo-relative permeability approach [8, 39, 6]. The pseudo-relative permeability approach computes nonlinear relative permeability functions. These nonlinear approaches are known to lack robustness and are process dependent [23, 24]. To overcome these difficulties, one needs a better understanding of nonlinear upscaling methods for space-time heterogeneous problems. Nonlinear upscaling methods for scale separation cases are rigorously treated in [43, 27]. Among these approaches, some deal with problems that have both space and time heterogeneities.

Our proposed approaches take their origin in the Constraint Energy Minimizing Generalized Multiscale Finite Element Method (GMsFEM) and Nonlocal Multi-Continua upscaling, which are related. The main idea of these approaches is to use multiple macroscopic parameters to represent the solution over each coarse-grid block. We refer to these degrees of freedom as continua, which are important for achieving a high order accuracy. We note that generalized continua concepts are also introduced in computational mechanics [28], which include generalized continuum theories (e.g., [28]), computational continua framework (e.g., [35]), and other approaches. Computational continua ([35, 29]), which use nonlocal quadrature to couple the coarse scale system stated on unions of some disjoint computational unit cells, are introduced for non-scale-separation heterogeneous media. In [34, 33, 30], the computational continua with model reduction technique is combined.

An important step that connects multiscale methods and upscaling techniques includes using basis functions such that the resulting degrees of freedom have physical meanings, typically averages of the solution. For nonlinear problems, using linear basis functions is not very suitable. The local problems are nonlinear problems. For this reason, in our first work [17], we provided a framework for NLMC for stationary problems. In this paper, we provide a unified framework for nonlinear NLMC for problems with space-time heterogeneities, analysis, and machine learning based simplified local solves.

In Figure 1, we illustrate the main steps of our approach. Below, we briefly describe them. In the first step, we identify continua in each coarse block. This is done with the help of test functions, which can separate the features that can not be localized within the region of influence (oversampling region designated with green color in Figure 1). For nonlinear problems, each continua is defined by a corresponding test function. Continua play the role of macroscale variables. In our examples, macroscale variables are average solution values in some selected heterogeneous regions (such as channels).

In Step 2, once we identify the continua, we use oversampling regions to define downscaling maps. The oversampling region represents the region of influence and thus, the macroscopic parameter interactions are defined within oversampling regions. The local nonlinear problems are formulated in the oversampled regions using constraints. However, these computations are expensive and require appropriate local problems. Instead, we propose to use local space-time models of the original PDEs and perform many tests with various boundary conditions and sources. These local solutions are used to train macroscopic parameters as a function of multiple macroscale continua variables. For machine learning, we use deep learning algorithms, which allow approximating complex multi-continua dependent functions.

In Step 3, we seek a coarse-grid solution (the values in each continua) such that the downscaled global fine-scale solution satisfies the variational formulation that uses the test functions defined in Step 1. An example of test functions that we use is piecewise constant functions in each subregions (defined as channels). Then, the macroscale variables are average solutions defined in these subregions. The corresponding downscaled maps represent the local fine-grid solutions given these constraints. The global coarse-grid formulation can be thought as a mass balance equation formulated for each continua.

The main contributions of this paper are the following:

• •

Novel upscaled model for space-time;

• •

Unified framework using test functions;

• •

Easy local problems and machine learning calculations;

• •

Numerical results that uses machine learning and nonlinear upscaled models.

In the paper, we present an analysis of our approach for a model problem, which consists of heterogeneous p-Laplacian ($p=2$). This model problem requires nonlinear upscaling and some oversampling in order to show an optimal convergence of our proposed approach.

In conclusion, the paper is organized as follows. In Section 2, we give some preliminary results of the nonlocal multicontinua approach. In Section 3, we present our approach, which uses the space-time nonlocal multicontinua approach. In this section, we present examples and convergence results. The numerical results are presented in Section 4.

2 Overview of NLMC methods

In this section, we will give a brief overview of the NLMC method for linear problems [16]. Our goal is to summarize the key ideas and motivate our new space-time nonlinear NLMC method. We consider a model elliptic equation with a heterogeneous coefficient

[TABLE]

Here $\kappa$ is the heterogeneous field, $f$ is a given source and $\Omega$ is the physical domain.

The NLMC method is defined on a coarse mesh, $\mathcal{T}^{H}$, of the domain $\Omega$. We write $\mathcal{T}^{H}=\bigcup\{K_{i}\;|\;i=1,\cdots,N\}$, where $K_{i}$ denotes the $i$-th coarse element and $N$ denotes the number of coarse elements in $\mathcal{T}^{H}$. For each coarse element $K_{i}$, we define an oversampled region $K_{i}^{+}$, which is obtained by enlarging the coarse block $K_{i}$ by a few coarse grid layers. We will also denote $K_{i}^{+}=K_{i,l}$ when the oversampling region is obtained by enlarging $K_{i}$ by $l$ coarse grid layers. See Figure 2 for an illustration of coarse grid and oversample region. In particular, a structured coarse grid is shown with boundaries of coarse elements are denoted red. A coarse cell $K_{i}$ is denoted green and its oversampled region $K_{i}^{+}$ obtained by enlarging $K_{i}$ by one coarse grid layer is enclosed by black lines.

The NLMC method consists of three main ingredients:

Choice of continua. 2. 2.

Local basis functions. 3. 3.

Global coupling.

For each coarse element $K_{i}$, we will identify multiple continua corresponding to various solution features. This can be done via a local spectral problem or a suitable weight function. Using the definition of continua, we will define a set of local basis functions by solving some local problems on oversample regions. Then, the final NLMC system is defined using these multiscale basis functions and a suitable variational formulation. In the following, we will discuss these concepts in detail.

Now we will specify the definition of continuum that is used in our studies. For each coarse block $K_{i}$, we will identify a set of continua which are represented by a set of auxiliary basis functions $\phi^{j}_{i}$, where $j$ denotes the $j$-th continuum. There are multiple ways to construct these functions $\phi^{j}_{i}$.

One way is to follow the idea proposed in CEM-GMsFEM [18]. In this framework, the auxiliary basis functions $\phi^{j}_{i}$ are obtained as the dominant eigenfunctions of a local spectral problem defined on $K_{i}$. These eigenfunctions can capture the heterogeneities and the contrast of the medium. We can also follow the framework in the original NLMC method [16], 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 $\phi^{j}_{i}$ are piecewise constant functions, namely, they equal one within one fracture network and zero otherwise. Moreover, one can define the continua by using properties of the heterogeneous media. In this case, the auxiliary basis functions are piecewise constant functions defined with respect to a partition of the coarse cell $K_{i}$, such as the medium coefficients have a bounded contrast in each subregion [47].

Once the auxiliary basis functions $\phi^{j}_{i}$ are specified, we can construct the required basis functions. The idea generalizes the original energy minimization framework in CEM-GMsFEM. First, we denote the space of auxiliary basis functions as $V_{aux}$. Consider a given coarse element $K_{i}$ and a given continuum $j$ within $K_{i}$. We will use the corresponding auxiliary basis function $\phi^{j}_{i}$ to construct our required multiscale basis function $\psi^{j}_{i}$ by solving a problem in an oversampled region $K_{i}^{+}$. Specifically, we find $\psi^{j}_{i}\in H^{1}_{0}(K_{i}^{+})$ and $\mu\in V_{aux}$ such that

[TABLE]

where $\delta_{im}$ denotes the standard delta function and $\tilde{\kappa}$ is a weight function. We remark the function $\mu$ serves as a Lagrange multiplier for the constraints in the second equation of (2). We also remark that the basis function $\psi^{j}_{i}$ has mean value one on the $j$-th continuum within $K_{i}$ and has mean value zero in all other continua in all coarse elements within $K_{i}^{+}$. In practice, the above system (2) is solved in $K_{i}^{+}$ using a fine mesh, which is typically a refinement of the coarse grid. See Figure 2 for an illustration.

Finally, we can derive the NLMC system. Let $V_{ms}$ be the space spanned by the basis functions $\{\psi^{j}_{i}\}$. We will represent the approximate solution $u_{ms}\in V_{ms}$ as a linear combination of basis functions, namely,

[TABLE]

Then, we will find $u_{ms}$ by the following variational formulation

[TABLE]

This variational formulation results in the following upscaled model for the solution $U=(U^{j}_{i})$:

[TABLE]

where the upscaled stiffness matrix $A_{T}$ is defined as

[TABLE]

and the upscaled source term $F$ is defined as

[TABLE]

We remark that the nonlocal connections of the continua are coupled by the matrix $A_{T}$. We also remark that the local computation in (2) results from a spatial decay property of the multiscale basis function, see [18, 19, 14] for the theoretical foundation.

The above NLMC idea can be extended to nonlinear elliptic problems, resulting in a nonlinear NLMC method (14)-(15). See Section 3.4 for the derivation and the convergence analysis.

3 Nonlinear non-local multicontinua model

In this section, we present the nonlinear non-local multicontinua (NLMC) method. We will first give some general concept of the methodology in Section 3.1. Then, in Section 3.2, we give some illustrative examples including linear problems and pseudomonotone problems. The main methodological details of the method are presented in Section 3.3. Finally, we present a convergence analysis of the method for a model elliptic problem in Section 3.4.

3.1 General concept

We will first present some general concepts of our nonlinear NLMC using the following model nonlinear problem

[TABLE]

where $G$ is a nonlinear operator that has a multiscale dependence with respect to space (and time, in general) and $M$ is a linear operator. In the above equation, $U$ is the solution and $g$ is a given source term. Our method has three key ingredients, namely, the choice of continua, the construction of local downscaling map and the construction of the coarse scale model. We will summarize these concepts in the following.

• •

The choice of continua

The continua serve as our macroscopic variables in each coarse element. Our approach uses a set of test functions to define the continua. To be more specific, we consider a coarse element $K_{i}$. We will choose a set of test functions $\{\psi_{i}^{(j)}(x,t)\}$ to define our continua, where $j$ denotes the $j$-th continuum. Using these test functions, we can define our macroscopic variables as

[TABLE]

where $\langle\langle\cdot,\cdot\rangle\rangle$ is a space-time inner product.

• •

The construction of local downscaling map

Our upscale model uses a local downscaling map to bring microscopic information to the coarse grid model. The proposed downscaling map is a function defined on an oversampling region subject to some constraints related to the macroscopic variables. In time-dependent problems, the oversampling region can be regarded as a zone of influence for coarse-grid variables defined on the target coarse block $K_{i}$. More precisely, we consider a coarse element $K_{i}$, and an oversampling region $K_{i}^{+}$ such that $K_{i}\subset K_{i}^{+}$. Then we find a function $\phi$ by solving the following local problem

[TABLE]

The above equation (5) is solved subjected to constraints defined by the following functionals

[TABLE]

This constraint fixes some averages of $\phi$ with respect to $\psi_{i}^{(j)}(x,t)$. We remark that the function $\mu$ serves as the Lagrange multiplier for the above constraints. This local solution builds a downscaling map

[TABLE]

• •

The construction of coarse scale model

We will construct the coarse scale model using the test functions $\{\psi_{i}^{(j)}(x,t)\}$ and the local downscaling map. Our upscaling solution $U^{ms}$ is defined as a combination of the local downscaling maps. To compute $U^{ms}$, we use the following variational formulation

[TABLE]

The above equation (6) is our coarse scale model.

We would like to briefly summarize above steps. The first step defines multicontinua, which play the role of macroscopic variables. They are critical in multiscale modeling and need to be defined apriori. The second step constructs downscaling maps and can be computationally intensive. We will propose a machine learning technique in combination with solving local problems of the original equation subject to various boundary conditions. From here, the macroscale fluxes will be defined as a function of macroscopic variables in oversampled regions. This high dimensional functions will be learned using machine learning techniques during coarse-grid solution step (Step 3). Next, we will give some examples (Section 3.2) and then present a more detailed description of the algorithm (Section 3.3).

3.2 Examples

We will present two model problems, and discuss how our nonlinear NLMC is applied.

3.2.1 Linear case

In this section, we will construct our upscaling model for a case that $G$ is a linear operator. We will follow the general concepts in Section 3.1. First, we discuss the choice of continua. For each coarse element $K_{i}$, we consider a set of test functions $\{\psi^{(j)}_{i}(x,t)\}$ defined for $x\in K_{i}$. Here the index $j$ denotes the $j$-th continuum. One choice of these test functions is a set of piecewise constant functions. Another choice of these test functions is the first $j$ dominant eigenfunctions of an appropriate spectral problem.

Next, we discuss the construction of the local downscaling map. We fix a continuum $\psi^{(j)}_{i}(x,t)$ in the coarse region $K_{i}$. Let $K_{i}^{+}$ be an oversampling region. With the assumption that $G$ is linear, we can represent the downscaling map, denoted by $\phi^{(j)}_{i}$, as a linear combination of some generic local solutions $\{\phi^{(j,l)}_{i,m}\}$. To find these functions $\{\phi^{(j,l)}_{i,m}\}$, we solve the following

[TABLE]

on the oversample region $K_{i}^{+}$, where $\langle\langle\cdot,\cdot\rangle\rangle$ is an inner product and $\mu^{(j,l)}_{i,m}$ plays the role of Lagrange multiplier. Using these functions $\{\phi^{(j,l)}_{i,m}\}$, we can represent the local downscaling map $\psi^{(j)}_{i}(x,t)$ as

[TABLE]

Since $G$ is linear, we have

[TABLE]

Let $\{\chi_{i}\}$ be a set of partition of unity functions corresponding to the partition $\{K_{i}^{+}\}$ of the domain $\Omega$. The final upscale solution is then defined as the combination $\phi:=\sum_{i}\sum_{j}\chi_{i}\phi_{i}^{(j)}$. Using the test functions $\psi_{i}^{(j)}$, we can compute the macroscopic value $\{U_{m}^{(l)}\}$ by the following variational formulation

[TABLE]

3.2.2 Pseudomonotone case

Next, we consider another example for which $G$ is a pseudo-monotone operator. In this case, to compute the downscaling map, $\mathcal{F}^{ms}$, we will need to solve the following local problem: find $\mathcal{F}^{ms}(U)$ and $\mu$ such that

[TABLE]

The coarse grid system is then defined as

[TABLE]

3.3 More details of general framework

In this section, we give the details of our nonlinear NLMC framework. We consider the following model problem of finding $u\in V$ such that

[TABLE]

with $u(\cdot,0)=0$, where $L$ is a nonlinear differential operator, $T>0$ is a fixed time and $V$ is a suitable function space. We use a different notation for nonlinear differential operator as in (4) to simplify the notations, and our methodology remains applicable to the problem described by (4).

Next, we discuss the mesh. We assume that $\Omega$ is partitioned by a coarse mesh $\mathcal{T}_{H}$ (see Figure 2) with mesh size $H>0$ and $(0,T]$ is partitioned into coarse time intervals denoted as $\mathcal{T}_{T}=\{(t_{i},t_{i+1}]\}$. A space-time element $K^{(n,i)}$ is then defined by $K_{i}\times(t_{n},t_{n+1}]$ for a coarse cell $K_{i}\in\mathcal{T}_{H}$ and the $n$-th time interval $(t_{n},t_{n+1}]$. The construction of our nonlinear NLMC method follows the three steps explained in Section 3.1.

Approximation by global basis

The discussion of our method starts with the use of global basis functions. In this case, the basis functions are global in space and in time. The motivation of this follows from the global basis of CEM-GMsFEM [18], for which coarse grid convergence is obtained.

• •

Choice of continua

The continua is defined using a set of test functions. Consider a space-time element $K^{(n,i)}$, we will introduce a set of test functions $V_{aux}=\{\psi_{j}^{(n,i)}\}$ which corresponding to different continua of the problem. We notice that $\psi_{j}^{(n,i)}$ is supported in $K^{(n,i)}$. We let $N_{c}$ be the number of such test functions. Then we will define macroscopic variables by

[TABLE]

where $s(\cdot,\cdot)$ is a weighted $L^{2}$ inner product with weighting function $\tilde{\kappa}$ such that $c_{0}H^{-1}\leq\tilde{\kappa}\leq c_{0}H^{-1}$. Note that this condition for the weighting function is motivated by the weighting function used in CEM-GMsFEM.

• •

Global downscaling map

We will define a downscaling map. This downscaling map will give a function defined globally in space and in time with constraints defined using a given set of macroscopic values. More precisely, we fix a set of macroscopic values $\{U_{j}^{(n,i)}\}$. We will then define a function $F=(F_{1},F_{2})$ such that $F_{1}\in V$ and $F_{2}\in V_{aux}$. These functions are obtained by solving

[TABLE]

We notice that the global function $F_{1}$ has macroscopic values equal to the given values $\{U_{j}^{(n,i)}\}$ and the function $F_{2}$ serves as the Lagrange multiplier for these constraints.

• •

Coarse grid model

Next, using the downscaling map, we can define the global coarse grid problem as: finding $U\in\mathbb{R}^{N_{c}}$ such that

[TABLE]

Then, the global numerical solution $u_{glo}$ is defined by $u_{glo}=F_{1}(U)$.

Nonlinear NLMC method

Now we will present the nonlinear NLMC method. The key ingredient is that we will replace the global downscaling map above by a local downscaling map.

• •

Local downscaling map

We will introduce the localized downscaling operator $F_{ms}=(F_{ms,1},F_{ms,2})$. Consider a space-time element $K^{(n,i)}$. We define a space-time oversampling region $K_{+}^{(n,i)}=K_{i}^{+}\times(t_{n}^{-},t_{n+1}]$ where $t_{n}^{-}<t_{n}$. We will then define a function $F_{loc}^{(n,i)}=(F_{loc,1}^{(n,i)},F_{loc,2}^{(n,i)})$ such that $F_{loc,1}^{(n,i)}\in V(K_{+}^{(n,i)})$ and $F_{loc,2}^{(n,i)}\in V_{aux}(K_{+}^{(n,i)})$, where $V(K_{+}^{(n,i)})$ and $V_{aux}(K_{+}^{(n,i)})$ are restrictions of $V$ and $V_{aux}$ on $K_{+}^{(n,i)}$ respectively. These functions are obtained by solving

[TABLE]

Finally the localized downscale operator is defined by $F_{ms,p}(U)=\sum_{n,l}\chi^{(n,i)}F_{loc,p}^{(n,i)}(U)$ where $p=1,2$ and $\chi^{(n,i)}$ is a partition of unity such that $\sum_{n,i}\chi^{(n,i)}\equiv 1.$

• •

Coarse grid model

The coarse grid problem is then defined as: finding $U\in\mathbb{R}^{N_{c}}$ such that

[TABLE]

and the nonlinear NLMC solution $u_{ms}$ is defined by $u_{ms}=F_{ms,1}(U)$.

3.4 Error sources and analysis

In this section, we present a concept of the analysis for the method. We will use a simple monotone elliptic equation to illustrate the main ideas. We consider the following problem: find $u$ such that

[TABLE]

where $\kappa(x,v)$ is a heterogeneous function. The weak formulation of the above equation can be written as: find $u\in V=H_{0}^{1}(\Omega)$ such that

[TABLE]

where, for any open subset $\omega\subset\Omega$ of the domain, the operator $A_{\omega}$ is defined by

[TABLE]

We will assume that the heterogeneous function $\kappa(x,v)$ satisfies the following two properties.

Assumption on $\kappa(x,v)$

If the vector field $v=0$, then $\kappa(x,v)=0$. 2. 2.

Lipschitz continuity with respect to $v$:

We assume there exist a function $\overline{\kappa}\in L^{\infty}(\Omega)$ such that

[TABLE] 3. 3.

Monotonicity:

We assume that the following coercivity condition holds

[TABLE]

Next, for any open subset $\omega\subset\Omega$ of the domain, we define two inner products $a_{\omega}(\cdot,\cdot)$ and $s_{\omega}(\cdot,\cdot)$ as follows

[TABLE]

where $\tilde{\kappa}(x)=\overline{\kappa}\sum_{i}|\nabla\chi_{i}|^{2}$ and $\{\chi_{i}\}_{i=1}^{N}$ is a set of partition of unity functions corresponding to the coarse mesh such that $0\leq\chi_{i}\leq 1$. The norms $\|\cdot\|_{a(\omega)}$ and $\|\cdot\|_{s(\omega)}$ corresponding to these inner products are defined as

[TABLE]

respectively. To simplify the notation, we use $A$, $\|\cdot\|_{a}$ and $\|\cdot\|_{s}$ to denote $A_{\Omega}$, $\|\cdot\|_{a(\Omega)}$ and $\|\cdot\|_{s(\Omega)}$ respectively.

In the following Lemma, we will show that the operator $A_{\omega}$ satisfies some coercivity and continuity properties.

Lemma 1**.**

For $\omega\subset\Omega$, $u,v,w\in H^{1}(\omega)$, we have

[TABLE]

and

[TABLE]

Moreover, we have

[TABLE]

Proof.

By the assumption (13), we have

[TABLE]

which proves the first inequality. By the assumption (12), we have

[TABLE]

which gives the second inequality. Recall that $\kappa(x,0)=0$. Thus we have

[TABLE]

which shows the third inequality. This completes the proof of this lemma. ∎

We next prove the following technical result.

Lemma 2**.**

For $\omega\subset\Omega$, $u,v\in H^{1}(\omega)$, we have

[TABLE]

Proof.

Notice that $\nabla(\chi_{i}v)=v\nabla\chi_{i}+\chi_{i}\nabla v$. Thus, we have

[TABLE]

We will first estimate the term $\int_{\omega}v\kappa(x,\nabla u)\cdot\nabla\chi_{i}$. By the Cauchy-Schwarz inequality, we have

[TABLE]

and, by using (12), we have

[TABLE]

Combining the above, we have

[TABLE]

To estimate the second term $\int_{\omega}\chi_{i}\kappa(x,\nabla u)\nabla v$, we use the fact that $|\chi_{i}|\leq 1$ and assumption (12) to obtain

[TABLE]

This completes the proof of this lemma. ∎

In the following, we will formulate our nonlinear NLMC method for the equation (11). The coarse scale degrees of freedom (continua) $\overline{U}$ of the solution is defined as

[TABLE]

for some $\mu_{i,j}\in L^{\infty}(\Omega)$ where $\mu_{i,j}|_{K_{m}}=0$ if $m\neq i$. The auxiliary space $V_{aux}$ is then defined as

[TABLE]

We remark that the functions in $V_{aux}$ defines the continua. In particular, $\mu_{i,j}$ defines the $j$-th continuum in the coarse cell $K_{i}$.

Next, to construct the numerical upscaling equation for (11), we will define a global downscaling operator $F_{1}$ such that $F_{1}(\overline{U})\in H_{0}^{1}(\Omega)$ and

[TABLE]

Next, we will define a projection operator $\Pi:V\rightarrow V_{aux}$ such that

[TABLE]

The global solution $U_{glo}\in V_{aux}$ is defined by

[TABLE]

and the global downscaled solution $u_{glo}$ is defined as $u_{glo}=F_{1}(U_{glo})$.

Approximation by global basis

We summarize the main steps:

Find $U_{glo}\in V_{aux}$

[TABLE] 2. 2.

Define

[TABLE]

Next, we will construct the nonlinear NLMC method. For each $K\in\mathcal{T}_{H}$, we will define a local downscaling operator $F^{loc,K}$ such that $F_{1}^{loc,K}(\bar{U})\in H_{0}^{1}(K^{+})$ and

[TABLE]

The multiscale solution $U_{ms}\in V_{aux}$ is defined by

[TABLE]

and the downscaled multiscale solution $u_{ms}$ is defined as $u_{ms}=F_{1}^{ms}(U_{ms}):=\sum_{K}\chi_{K}F_{1}^{loc,K}(U_{ms})$ where $\chi_{K}$ is a partition of unity such that $\sum_{K}\chi_{K}\equiv 1$ with $\sum_{K}|\nabla\chi_{K}|^{2}\leq C\sum_{i}|\nabla\chi_{i}|^{2}$ and $\text{supp}\{\chi_{K}\}\subset K^{+}=\cup_{\overline{K}_{i}\cap\overline{K}\neq\emptyset}\overline{K}_{i}.$

Nonlinear NLMC method

We summarize the main steps:

Find $U_{ms}\in V_{aux}$

[TABLE] 2. 2.

Define

[TABLE]

The analysis of our scheme is based on three assumptions. We summarize them below.

Assumption 1: For all $K\in\mathcal{T}_{H}$, $v\in V(K)$, we have

[TABLE]

Assumption 2: For $K\in\mathcal{T}_{H}$, $v_{aux}\in V_{aux}(K)$, there exist a function $w\in H_{0}^{1}(K)$ such that

[TABLE]

Assumption 3: There exist a $C_{\kappa}>0$ such that

[TABLE]

We will prove the following lemma for the stability of the downscale map.

Lemma 3**.**

By assumption 2, we have

[TABLE]

and

[TABLE]

Proof.

First, by Lemma 1 and (14), we have

[TABLE]

and by (15), we have

[TABLE]

Therefore, we have

[TABLE]

By Assumption 2, there exist a function $w\in H_{0}^{1}(K)$ such that

[TABLE]

Hence, we have

[TABLE]

This shows the first required inequality. Using a similar argument, we can prove that

[TABLE]

This completes the proof of the lemma. ∎

In the following lemma, we will give an error bound for the solution $F_{1}(U_{glo})$.

Lemma 4**.**

Let $u$ be the solution of (11) and $F_{1}(U_{glo})$ be the solution of (16)-(17). We have

[TABLE]

Proof.

First of all, we note that $F_{2}(U_{glo})=\Pi(f\tilde{\kappa}^{-1})$. So, we have

[TABLE]

and

[TABLE]

Therefore, by (13), we have

[TABLE]

where the last inequality follows from Assumption 1. This completes the proof. ∎

In the next lemma, we give a localization result. To do so, we need some notations for the oversampling domain and the cutoff function with respect to these oversampling domains. For each coarse cell $K$, we denote $K_{m}^{+}\subset\Omega$ as the oversampling coarse region by enlarging $K$ by $m$ coarse grid layers. For $M>m$, we define $\chi_{M,m}\in\text{span}\{\chi_{i}\}$ such that $0\leq\chi_{M,m}\leq 1$ and

[TABLE]

Note that, we have $K_{m}^{+}\subset K_{M}^{+}$.

Lemma 5**.**

Assume $K_{M}^{+}$ is an oversampling region obtained by enlarging the coarse cell $K$ by $M$ coarse grid layers. Let $\eta_{i}=F_{i}(\overline{U})-F_{i}^{loc,K}(\overline{U})$. We have

[TABLE]

and

[TABLE]

Proof.

The first step of the proof is to show the following inequality

[TABLE]

To do so, we denote $f_{i}=F_{i}(\overline{U})$ and $g_{i}=F_{i}^{loc,K}(\overline{U})$. By (13), we obtain

[TABLE]

Recalling that $\eta_{i}=f_{i}-g_{i}$. We notice that

[TABLE]

Therefore, we have

[TABLE]

Next, we define $K_{m}^{{}^{\prime}}=K_{m+1}^{+}\backslash K_{m}^{+}$ and obtain

[TABLE]

Therefore, we have

[TABLE]

Next we will estimate $\|\eta_{2}\|_{s}$. By Assumption 2, for $K\in\mathcal{T}_{H}$, there exist a $v\in H_{0}^{1}(K)$ such that

[TABLE]

Thus, for $K\subset K_{M}^{+}$, we have

[TABLE]

Combining (23) and (24), we have

[TABLE]

This shows (22).

By using (22), we have

[TABLE]

and we therefore obtain

[TABLE]

This gives the first required inequality. The second required inequality follows from (24). This completes the proof of tis lemma. ∎

The following result gives an error estimate for our nonlinear NLMC solution.

Theorem 1**.**

Consider the oversampling domain $K_{M}^{+}$ obtained by enlarging $K$ by $M$ coarse cell layers. Let $u$ be the solution of (11) and $F^{ms}_{1}(U_{ms})$ be the solution of (18)-(19). Then we have

[TABLE]

where

[TABLE]

Moreover, if M\sim O\Big{(}\log(H^{-1})+\log(C_{\kappa})\Big{)} and $H\leq\cfrac{1}{2}$, then we have

[TABLE]

Proof.

We will analyze the error by first separating the error into three parts as follows

[TABLE]

By Lemma 4, we have

[TABLE]

and by Lemma 5, we have

[TABLE]

By Lemma 3, we have

[TABLE]

and

[TABLE]

Therefore, we obtain

[TABLE]

Next, we will estimate the term $\|F_{1}(U_{glo})-F_{1}(U_{ms})\|_{a}$. By Lemma 1 and Assumption 3, we have

[TABLE]

and by Lemma 5, we have

[TABLE]

Therefore, we have

[TABLE]

Combining the above results, we obtain

[TABLE]

To show the second part of the theorem, we notice that

[TABLE]

If $M$ is large enough such that

[TABLE]

then we have

[TABLE]

and

[TABLE]

This completes the proof of the theorem.

∎

4 Numerical results

In this section, we present numerical results for the proposed method. In our examples, we will use simplified local problems to compute macroscale parameters. These local computations will involve machine learning algorithms. We consider following model problems in fractured and heterogeneous porous media:

• Test 1: Nonlinear flow problem (unsaturated flow problem)
• Test 2: Nonlinear transport and flow problem (two-phase flow problem)

We solve model problem in $\Omega=[0,1]\times[0,1]$ with no flux boundary conditions. Heterogeneous porous matrix permeability and location of the source terms and fracture position are depicted in Figure 3. We set source terms $f^{\pm}=\pm q$, $q=10^{3}$. We use $10\times 10$ coarse grid and $160\times 160$ fine grid.

*Test 1. * We consider the solution of the nonlinear equation in fractured heterogeneous media. For the nonlinear coefficients, we use $k^{\alpha\beta}(x,u)=k_{s}(x)k_{r}(u)$ with $k_{r}(u)=\exp(-a|u|)$, $a=0.1$ ($\alpha,\beta=m,f$). We set $c^{m}=1$, $c^{f}=0$, $k_{s}^{f}=10^{6}$ and $T_{max}=10^{-3}$ with 20 time steps.

*Test 2. * We consider the solution of the two-phase flow problem in fractured and heterogeneous porous media. For nonlinear coefficients, we set $\lambda^{w}(s)=s^{2}$ and $\lambda^{n}(s)=(1-s)^{2}$. We set $\phi^{\alpha}=1$ ($\alpha=m,f$), $k^{f}=10^{3}$ and $T_{max}=6.3\cdot 10^{-5}$ with 700 time steps.

Each sample $X_{l}$ contains the information about heterogeneous permeability and fracture positions up to the fine grid resolution in local domain, coarse grid mean value of the solution in oversampled local domain

[TABLE]

and output

[TABLE]

Each dataset is divided into training and validation sets with $80:20$ ratio.

For the training of the neural networks, we use a global dataset, where we extract local information from the fine grid calculations on the global domain $\Omega$. We train four neural networks for each type of transmissibility: $NN_{1}$ for horizontal coarse edges for matrix-matrix flow, $NN_{2}$ for vertical coarse edges s for matrix-matrix flow, $NN_{3}$ for matrix - fracture flow and $NN_{4}$ for fracture - fracture flow. For calculations, we use $150$ epochs with a batch size $N_{b}=90$ and Adam optimizer with learning rate $\epsilon=0.001$. For accelerating of the training process of the multi-input CNN, we use GPU. We use $3\times 3$ convolutions and $2\times 2$ maxpooling layers with RELU activation for $X^{k}$ and $X^{f}$, and $3\times 3$ convolutions with RELU activation for $X^{\overline{p}^{m}}$. For each input data, we have 2 layers of CNN with one final fully connected layer. Convolution layer contains 8 and 16 feature maps for $X^{k}$ and $X^{f}$; and 4 and 8 feature maps for $X^{\overline{p}^{m}}$. We use dropout with rate 10 % in each layer in order to prevent over-fitting. Finally, we combine CNN output and perform two additional fully connected layers with size 200 and 1(one final output). Presented algorithm is used to learn dependence between multi-input data and upscaled nonlinear transmissibilities.

For error calculation on the dataset, we used mean square errors, relative mean absolute and relative root mean square errors

[TABLE]

where $Y_{i}$ and $\tilde{Y}_{i}$ denotes reference and predicted values for sample $X_{i}$ Learning performance for neural networks are presented in Table 1 for Test 1 and Test 2. We observe a good convergence with small error for each neural network.

Next, we consider errors between solution of the coarse grid problem with the reference and predicted upscaled transmissibilities. To measure difference between reference solution and coarse grid solution, we compute relative $L_{2}$ error

[TABLE]

where $u=p,s$, $\overline{u}^{fine}$ is the reference solution (mean value on coarse grid of the fine grid solution) and $\overline{u}$ is the solution on the coarse grid. In Figure 4, we depict solution of the problem for Test 1 on the fine grid, coarse grid upscaled solution using classic approach and for new method presented ($u^{fine}$, $\overline{u}^{fine}$, $\overline{u}^{UP}$ and $\overline{u}^{NL}$). We have $e(\overline{u}^{UP})=11.773\%$ and $e(\overline{u}^{NL})=2.155\%$ at final time.

In Figure 5, we depict the solution of the problem for Test 2. On the first column, we depict a reference fine grid solution ($s^{fine}$, $p^{fine}$), mean value on coarse grid of the fine grid solution ($\overline{s}^{fine}$, $\overline{p}^{fine}$) on the second column, coarse grid solution using upscaling method ($\overline{s}^{UP}$, $\overline{p}^{UP}$) on the third column and coarse grid solution using nonlinear nonlocal machine learning method ($\overline{s}^{NL}$, $\overline{p}^{NL}$) on the fourth column. On the first, second and third rows, we show a saturation for time $t_{m}$, $m=300,700$ and on fourth row, we have pressure for time $t_{m}$, $m=700$. Fine grid (reference) solution is performed using finite volume approximation with embedded discrete fracture model, where for error calculations, we used a mean values of the reference solution on the coarse grid, $\overline{p}^{fine}$ and $\overline{s}^{fine}$. On the last column of the Figure 5, we depict a coarse grid solution using nonlinear nonlocal transmissibilities that calculate based on the machine learning approach. For machine learning approach, we have $e(\overline{p}^{NL})=0.281\%$, $e(\overline{s}^{NL})=3.512\%$, and for upscaling $e(\overline{p}^{UP})=14.063\%$, $e(\overline{s}^{UP})=13.354\%$ at final time $t_{m}$, $m=700$.

5 Conclusions

In the paper, we present a general nonlinear upscaling framework for nonlinear differential equations with multiscale coefficients. The framework is built on nonlinear nonlocal multi-continuum upscaling concept. The approach first identifies test functions for each coarse block, which are used to identify macroscale variables (called continua). In the second stage, we solve nonlinear local problems in oversampled regions with some constraints defined via test functions. Simplified local problems are proposed for numerical results. Deep learning algorithms are used to approximate the nonlinear fluxes that are derived in nonlinear upscaling. In the final stage, macroscale formulation is given and it seeks the values of macroscopic variables such that the downscaled field solves the global problem in a weak sense defined using the test function. We present an analysis of our approach for an example nonlinear problem. We present numerical results for several porous media applications, including two-phase flow and transport.

Acknowledgements

The research of Eric Chung is partially supported by the Hong Kong RGC General Research Fund (Project numbers 14304217 and 14302018) and CUHK Faculty of Science Direct Grant 2018-19.