A nonlinear multigrid method for inverse problem in the multiphase porous media flow

Applied Mathematics and Computation 320 (2018) 271–281

Contents lists available at ScienceDirect

Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc

A nonlinear multigrid method for inverse problem in the multiphase porous media flow Tao Liu School of Mathematics and Statistics, Northeast University at Qinhuangdao, Qinhuangdao 066004, PR China

a r t i c l e

i n f o

Keywords: Inverse problem Nonlinear multigrid Multiphase porous media flow

a b s t r a c t In this paper, we consider a parameter identification problem for the nonlinear convection–diffusion equation in the multiphase porous media flow. A nonlinear multigrid method is proposed for the recovery of permeability. This method works by dynamically adjusting the objective functionals at different grids so that they are consistent with each other, and ultimately reduce, the finest grid objective functional. In this manner, the nonlinear multigrid method can efficiently compute the solution to a desired fine grid inverse problem. Numerical results illustrate that the proposed multigrid approach both dramatically reduces the required computation and improves the reconstructed image quality. © 2017 Elsevier Inc. All rights reserved.

1. Introduction The multiphase flow in porous media has received considerable attention in recent years because of increasing applications in science and engineering. Applications of porous media flow include Biology [1], Mechanics [2], Geosciences [3], Civil Engineering [4]. In many engineering problems, such as those found in reservoir simulation, the oil and gas flow within the reservoir is often modeled as porous media flow. In [5], Espedal and Karlsen considered the saturation equation in the fractional flow formulation of the two-phase porous media flow equations, which is actually a convection dominated, degenerate convection–diffusion equation. Since then, various numerical techniques have appeared in the literatures of the inverse problem in porous media flow. One of the most utilized techniques for this problem is the ensemble Kalman Filter [6], and recently Hoel et al. [7] developed a multilevel version of this method. In [8], a contamination source identification problem in constant porous media flow is solved with a hierarchical Bayesian computation method. In [9], nuclear magnetic resonance imaging data are used in an inverse problem methodology to estimate the porous media flow functions. Berre et al. [10] solved the multi-level parameter structure identification problem for two-phase porous media flow problems using flexible representations, and Krüger et al. [11] estimated a spatially dependent diffusion function in porous media flow utilizing pressure and fluid rate data. Liu [12] and Nilssen et al. [13] respectively applied the wavelet multiscale-homotopy method and the augmented Lagrangian method to the reconstruction of a permeability field in the multiphase porous media flow. The inverse problem discussed in this paper is to identify q(x) in the following nonlinear convection-diffusion equation in the multiphase porous media flow:

ut + ∇ · (ϕ , ψ ) − ∇ · (q(x ) · N (u )∇ u ) = s(x, t ),

E-mail address: [email protected] https://doi.org/10.1016/j.amc.2017.09.039 0 096-30 03/© 2017 Elsevier Inc. All rights reserved.

in  × (0, T ),

(1.1)

272

T. Liu / Applied Mathematics and Computation 320 (2018) 271–281

with the initial-boundary condition



u(x, 0 ) = φ (x ), in , u(x, t ) = 0, on ∂  × (0, T ),

(1.2)

where u is the concentration field, q is the permeability coefficient, x is the space variable, t is the time variable, N is the positive nonlinear function in the diffusion term, s is the piecewise smooth source function, ∇ · (ϕ , ψ ) is defined by ∇ · (ϕ , ψ ) = ∂∂x ϕ (u ) + ∂∂y ψ (u ), and ϕ , ψ are the nonlinear S-shaped flux functions of Buckley–Leverett type in the x, y directions respectively. Eq. (1.1) is related to the multiphase flow in porous media. The immiscible displacement of oil by water in a porous media with zero gravity effects can be described by a set of nonlinearly coupled partial differential equations [5]

 ∇ · υ = s1 ( x, t ), υ = −q(x )κ (x, u )(∇μ − θ (u )∇ζ ), η (x )ut + ∇ · (φ (u )υ + ϕ (u )q(x )∇ζ ) − ∇ · (q(x )N (u )∇ u ) = s2 (x, t ),

(1.3)

where s1 and s2 are the injection and production wells, q is the permeability, υ is the total Darcy velocity, κ is the total mobility of the phases, μ is the global pressure, θ is the density of the wetting phase, ζ is the height, η is the porosity, N is the nonlinear diffusion function, φ is the nonlinear S-shaped fractional flow function, ϕ is defined by ϕ (u ) = (θ (u ) − ρ (u ))φ (u ) (x, u ), and ρ , are respectively the density and phase mobility of the nonwetting phase. Eq. (1.1) is similar to Eq. (1.3) except for the convection term and the time derivative term. The time derivative terms in these two equations are equal if we assume that η(x) ≡ 1. The difference in the convection terms is that Eq. (1.1) has no permeability dependence and no varying coefficient. The inverse problem discussed in this paper can be viewed as a parametric data-fitting problem, and we can formalize such a problem in the framework of optimization where a functional defined in terms of discrepancy between observed and computed data is minimized over a model space. Generally, such problem is very challenging to solve, since: (1) Ill-posedness: The solution does not depend continuously on the observed data. A minor disturbance of the observed data may cause large change on the solution of the inverse problem. (2) Nonlinearity: Nonlinear dependence of the observed data with respect to the permeability field to be identified causes the presence of numerous local minima. (3) Large computational cost: An inversion process needs tremendous forward computations, the computational cost is often very large. Therefore, the key problem is how to quickly find a stable solution. To overcome the three main difficulties, we utilize a regularized nonlinear multigrid method. Multigrid method is a specific form of multiresolution method that can be used to reduce the computational requirements and remove smooth error components of large numerical problems of partial differential equations [14–17]. This method works by recursively moving between different resolutions, thereby propagating information between coarse and fine scales. Multigrid method has been primarily used for solving forward problems, but more recently it has also been applied to a variety of inverse problems. Ascher and Haber [18] proposed an efficient multigrid method for the recovery of a coefficient function of an elliptic differential equation in three dimensions. Spiliopoulos et al. [19] applied a multigrid approach for the estimation of geometric anisotropy parameters from scattered spatial data that are obtained from environmental surveillance networks. Mewes et al. [20] studied the multigrid Monte Carlo solution of anisotropic seismic inversion. McCormick and Wade [21] applied multigrid methods to a linearized electrical impedance tomography problem, and Borcea [22] used a nonlinear multigrid approach to electrical impedance tomography based on a direct nonlinear formulation analogous to the full approximation scheme in nonlinear multigrid partial differential equation solvers. Ye et al. [23] formulated the multigrid approach directly in an optimization framework, and used the method to solve optical diffusion tomography problems. In related work, Oh et al. [24] formulated multigrid algorithms for the solution of a broad class of optimization problems arising from inverse problems. Importantly, both the approaches of Ye and Oh are based on the matching of objective functional derivatives at different grids. In this paper, we discuss the use of the nonlinear multigrid method for the solution of inverse problem for the nonlinear convection-diffusion equation in the multiphase porous media flow. The forward problem is discretized by the finitedifference method at different grids, and the inverse problem is formulated as a sequence of optimization problems. Then, we dynamically adjust the objective functionals in these optimization problems, in order to make them consistent, and ultimately reduce the objective functional at the finest grid. The motivation for using the multigrid method is that it can provide several unique advantages. First, some form of iteration is usually imperative in nonlinear problems, so multigrid methods are well adapted for solving nonlinear problems due to their iterative nature [25]. Second, since the coarse grid problems involve reduced number of variables, the number of local minima is enormously reduced and those that remain are further apart from each other. The reason lies in that by multigrid decomposition, lots of local minima are degenerated to trivial points, and the risk of trapping in a local minimum is reduced, thus multigrid methods can improve the inversion result. Third, multigrid methods can greatly reduce computation because both the forward and inverse problems are more coarsely discretized at coarser grids.

T. Liu / Applied Mathematics and Computation 320 (2018) 271–281

273

Outline of the paper is as follows, in Section 2 we present the numerical scheme that is used to discretize the forward problem, and then formulate the inverse problem in a discrete setting. In Sections 3 and 4, we respectively introduce the components of multigrid and design the concrete nonlinear multigrid method. Section 5 presents the results of numerical experiments. Finally, we conclude the paper with a summary of the work in Section 6. 2. Discretization and the inverse problem By using the finite-difference method, Eq. (1.1) can be discretized as

utk + ∇k · (ϕ , ψ ) − ∇k · (q(x ) · N (u )∇k u ) = sk (x, t ),

(2.1)

where the subscript k is the discretization parameter. For simplicity, we assume that  is the unit square, and use the following notation:

xi = ihx ,

i = 0, 1, . . . , m1 ;

tn = nht ,

yi = jhy , uni, j

n = 0, 1, . . . , m3 ;

j = 0, 1, . . . , m2 ;

= u ( xi , y j , tn );

sni, j = s(xi , y j , tn ),

where hx = m1 , hy = m1 , ht = mT are respectively spatial and time step lengths. 1 2 3 The forward and backward difference approximations of first derivative in the x direction are respectively defined by

Dx+ uni, j =

uni+1, j − uni, j hx

,

Dx− uni, j =

uni, j − uni−1, j hx

,

and the corresponding notations for the discrete derivatives in the y and t directions are used. The discretizations of the terms in Eq. (2.1) will be described in the following. The time derivative term is discretized as

utk = Dt− uni, j . By using the Engquist–Osher upwind scheme [26], the convection term is discretized as −1 ∇k · (ϕ (uni,−1 ), ψ (uni,−1 )) = Dx− ϕ EO (uni,−1 , uni+1 ) + Dy− ψ EO (uni,−1 , uni,−1 ), j j j ,j j j+1

where

1 2

−1 −1 ϕ EO (uni,−1 , uni+1 ) = (ϕ (uni,−1 ) + ϕ (uni+1 )) − j ,j j ,j

1 2

1 2

ψ EO (uni,−1 , uni,−1 ) = (ψ (uni,−1 ) + ψ (uni,−1 )) − j j+1 j j+1



−1 uni+1 ,j

uni,−1 j

1 2



|ϕ  ( ξ )|d ξ ,

uni,−1 j+1

uni,−1 j

|ψ  ( ξ )|d ξ .

The nonlinear diffusion term is discretized as x

y

∇k · (qi, j Ni,n j ∇k uni, j ) = Dx− (qyi+ 1 , j (N )ni+ 1 , j Dx+ uni, j ) + Dy− (qxi, j+ 12 (N )ni, j+ 1 Dy+ uni, j ), 2

2

2

where

1 1 (q 1 1 + qi− 12 , j+ 12 ), qyi+ 1 , j = (qi+ 12 , j+ 12 + qi+ 12 , j− 12 ), 2 i+ 2 , j+ 2 2 2 1 1 x n y n n n (N )i+ 1 , j = (N (ui+1, j ) + N (ui, j )), (N )i, j+ 1 = (N (uni, j+1 ) + N (uni, j )), 2 2 2 2 x

qi, j+ 1 = 2

and the discretized permeability q is defined by qi+ 1 , j+ 1 = q(xi+ 1 , y j+ 1 ). Eq. (2.1) can now be written as

2

2

2

Dt− uni, j + ∇k · (ϕ (uni,−1 ), ψ (uni,−1 )) − ∇k · (qi, j Ni,n j ∇k uni, j ) = sni, j , j j n = 1, 2, . . . , m3 ,

2

(2.2)

where u0i, j = φ (xi , y j ) and uni, j = 0 for (xi , yj ) ∈ ∂ . Using Eq. (2.2), we determine the density field u(x, y, t) with the known permeability q(x, y), such problem is denoted as a forward problem for the nonlinear convection–diffusion equation in the multiphase porous media flow. While the objective of the inverse problem is to identify the permeability q(x, y) from some observed data ud (xm , ym , t), m = 1, 2, . . . , M; t ∈ (0, T) of u(x, y, t). Generally, the inverse problem in the multiphase porous media flow can be viewed as a parametric data-fitting problem. It is possible to formulate such problem as an optimization problem where a functional defined by the discrepancy between observed and computed data is minimized over a model space. Now, to obtain this optimization problem, we first define a nonlinear operator equation by Eq. (2.2):

F ( q ) = u,

(2.3)

274

T. Liu / Applied Mathematics and Computation 320 (2018) 271–281

where q and u are vectors with the following form: 

q = (q 1 , 1 , q 1 , 3 , . . . , q 1 ,m2 − 1 , q 3 , 1 , p 3 , 3 , . . . , q 3 ,m2 − 1 , . . . , qm1 − 1 , 1 , qm1 − 1 , 3 , . . . , qm1 − 1 ,m2 − 1 ) , 2

2

2

2

2

2

2

2

2

2

2

2

2

2

2

2

2

2



3 3 3 u = (u1x1 ,y1 , u1x2 ,y2 , . . . , u1xM ,yM , u2x1 ,y1 , u2x2 ,y2 , . . . , u2xM ,yM , . . . , um , um , . . . , um ) . x1 ,y1 x2 ,y2 xM ,yM

Here the vector q is defined at half-points, because the permeability q(x) is defined by a single value within each reservoir simulator grid cell. Then, we can identify the permeability q from the observed data ud containing noise, by minimizing an output least-squares functional

min F (q ) − ud 2 ,

(2.4)

q

where  ·  is the Euclid norm. Actually, Eqs. (2.3) and (2.4) are just, respectively, the forward and inverse problems. The inverse problem is generally ill-posed, i.e., the parameter to be identified does not depend continuously on the data. So Tikhonov regularization [27] has to be used in order to compute a stable approximation of the minimizer, i.e., to minimize the associated least-squares problem

min K (q ) = F (q ) − ud 2 + μq2 , q

(2.5)

where μq2 is a regularization term, incorporated to impose stability or a priori information or both, and μ is a regularization parameter. In fact, the regularization term is constructed to provide the regularized solution. One intuitive explanation is that any nonzero feature that appears in the regularized solution increases the norm of q. Such features appear in the solution because they are necessary to fit the data. However, if noise is included in the data, and there is no point in fitting such noise exactly, obviously there might be many solutions adequately fit the data in the way that F (q ) − ud  ≤ δ . In general, all solutions with F (q ) − ud  ≤ δ are considered and among these solutions the one which minimizes q is selected. The reason for this is that the minimization of q can ensure that unneeded features will not appear in the regularized solution. 3. Components of multigrid In order to use the nonlinear multigrid method to identify the permeability q, we first require a sequence of discretization grids

{g }, g = 0, 1, . . . , G, where 0 is the finest discretization grid, which satisfies the resolution requirement for the permeability inversion in the multiphase porous media flow, and its discrete step length is (hx , hy ). The next coarser grid 1 should be defined by doubling the step length of 0 : (2hx , 2hy ). Then, the grid 2 can also be defined by doubling the step length of 1 : (22 hx , 22 hy ). In this way, all the grids g (g > 0) can be defined using a similar procedure, and its discrete step length is (2g hx , 2g hy ). The numerical method Eq. (2.2) resolving the forward problem should be stable at the coarsest grid G . For this stability issue, there are many good results in [26]. If we discretize the forward model at grid g , then we can obtain the nonlinear operator equation by Eq. (2.2):

F (g ) ( q (g ) ) = u (g ) ,

(3.1)

and identify the permeability q(g) by minimizing the following output least-squares functional

min K (g) (q(g) ) = F (g) (q(g) ) − ud(g) 2 + μ(g) q(g) 2 . q (g )

(3.2)

It is noteworthy that the forward model F(g) (q(g) ) and the optimizing functional K(g) (q(g) ) are both evaluated at grid g . This is very important because evaluation of the forward model at a coarse grid virtually reduces computation due to the reduction in the variable dimension, and evaluation of the optimizing functional at a coarse grid substantially improves the reconstructed image quality due to the reduction of local minima. The nonlinear multigrid method requires three components: restriction operator, prolongation operator, and relaxation operator. The restriction and prolongation operators are used to propagate the permeability information between coarse and fine grids, and the relaxation operator is used to update the permeability at a fixed grid. In this paper, these three operators are chosen as follows: g+1

(1) The restriction operators Ig



Ig+1 = g

1 1 4 1

1 1

g+1

: g → g+1 , g = 0, 1, . . . , G − 1, are chosen to be four point average:

. g

g

(2) The prolongation operators Ig+1 : g+1 → g , g = 0, 1, . . . , G − 1, are chosen to be constant prolongation. Constant prolongation means that an entry of q(g) at the fine discretization grid g is obtained by taking the value from the closest point of the coarse discretization grid g+1 .

T. Liu / Applied Mathematics and Computation 320 (2018) 271–281

275

(3) The relaxation operators Rg (q(g) , K (g) ), g = 0, 1, . . . , G, where q(g) is the initial value, and K(g) is the objective functional. A wide variety of numerical methods can be used as the relaxation operators. For example, common methods such as regularized Gauss–Newton method [28,29]



qk+1 = qk − F  (qk )∗ F  (qk ) + μI

−1 



· F  ( q k ) ∗ ( F ( q k ) − u d ) + μq k ,

and/or Levenberg Marquardt method [30]



qk+1 = qk − F  (qk )∗ F  (qk ) + μI

−1 



· F  ( qk )∗ ( F ( qk ) − ud ) ,

k = 0, 1, 2, . . .

k = 0, 1, 2, . . .

(3.3)

(3.4)

can be employed for solving the inverse problem at each grid. This makes our method particularly well suited to the solution of inverse problems. 4. Nonlinear multigrid method In the nonlinear multigrid method, the objective functional in the optimization problem is adjusted to be solved at each grid by using the solutions at both finer and coarser grids. To do this. we formulate a consistent set of coarse grid objective functionals to ultimately reduce the finest grid one. Let q(0) , q(g) respectively denote the finest grid permeability and a coarse grid representation of q(0) , and a coarser grid g+1 permeability q(g+1 ) can be obtained from q(g) by the relation q(g+1 ) = Ig q(g) . We adjust the objective functional K(g) ( · ) by appending an additional linear correction term as follows:

(g) (q(g) ) = F (g) (q(g) ) − u(g) 2 + μ(g) q(g) 2 − r(g) q(g) , K d

(4.1)

(0 ) (· ) = K (0 ) (· ), where is a row vector used to adjust the gradient of objective functional. At the finest grid, we need K ( 0 ) so r = 0. The next job is to explain how the objective functionals at each grid can be matched to produce a consistent (g ) solution. Specifically speaking, the recursive expressions of ud , μ(g) and r(g) are deduced such that the objective functionals at fine and coarse grids can match each other. Let q(g) be the current permeability approximation at grid g . To improve this approximation, we first perform an iterag+1 tion of the relaxation operator with the initial value Ig q(g) at the coarser grid g+1 : r(g)

(g) (g+1 ) q(g+1) ←− Rg+1 (Ig+1 ,K ), g q

(4.2)

where q(g+1 ) is the iteration result, and then use this result to correct the fine grid approximation by prolongating the change in the coarser grid permeability as (g ) q(g) = q(g) + Igg+1 (q(g+1) − Ig+1 ). g q

(4.3)

It is hoped that the new approximation q(g) should be at least as good as the old approximation q(g) , namely,

(g) (q(g) ). (g) (q(g) ) ≤ K K

(4.4)

But this will not be tenable unless the objective functionals are consistent. As a matter of fact, if we arbitrarily choose a set of objective functionals, the coarse grid correction will easily move the corrected result away from the optimal value. To avoid the problem of inconsistent objective functionals, we can use the condition that the fine and coarse grid objective functionals are equal within an additive constant for all values of q(g+1 )

(g+1) (q(g+1) ) ∼ (g) (q(g) + Ig (q(g+1) − Ig+1 q(g) )) + constant. K =K g g+1

(4.5)

(g+1 )

Then, our objective is to choose properly ud , μ(g+1 ) and r(g+1 ) , such that the coarse grid objective functional can match the fine grid objective functional as described in Eq. (4.5). First, we enforce the condition that the initial deviations between the forward model and observed data be the same at the coarse and fine grids, that is, (g ) F (g+1) (Ig+1 ) − ud(g+1) = F (g) (q(g) ) − ud(g) . g q

Hence, we can obtain



(4.6)



(g ) ud(g+1) = ud(g) − F (g) (q(g) ) − F (g+1) (Ig+1 ), g q

(4.7)

where the square bracket term compensates for the forward model mismatch between grids. Second, we use the condition that the initial regularization terms should be equal at the coarse and fine grids, giving (g ) 2 μ(g+1) Ig+1  = μ ( g ) q ( g ) 2 . g q

(4.8)

This yields the expression for μ(g+1 )

μ(g+1) =

q ( g ) 2 μ (g ) . ( g ) 2 Ig+1 g q

(4.9)

276

T. Liu / Applied Mathematics and Computation 320 (2018) 271–281

Fig. 1. The function of the term r(g+1) q(g+1) .

Finally, we introduce the necessary condition of convergence of the nonlinear multigrid inversion method [23,24], which is that the gradients of the coarse and fine grid objective functionals must be equal at the current parameter values of g+1 Ig q(g) and q(g) . More precisely, we will assume that the following condition is satisfied: (g ) ∇ K (g+1) (Ig+1 ) = ∇ K (g) (q(g) )Igg+1 , g q

(4.10)

g Ig+1

where the prolongation operator actually plays the role of the restriction operator because it multiplies the gradient vector on the right, and the role of the restriction operator (also known as decimation operator) is to propagate the information from the fine grid to the coarse grid. Then from Eq. (4.10), we can get the expression for r(g+1 ) as follows: (g ) r(g+1) = ∇ K (g+1) (Ig+1 ) − ∇ K (g) (q(g) )Igg+1 . g q

(4.11)

The above conditions (4.6), (4.8) and (4.10) are necessary to guarantee that the optimal value is a fixed point of the nonlinear multigrid inversion method. For their more detailed meaning, we refer the interested readers to the theorem in reference [24], which has given a set of sufficient conditions for monotone convergence of the multigrid inversion method. Fig. 1 graphically illustrates the necessity for using these conditions. Fig. 1a shows the case in which the gradients of the coarse and fine grid objective functionals are different at the current permeability approximation. In this case, the coarse grid objective functional cannot upper bound the value of the fine grid objective functional, and the iterated result may actually increase the value of the fine grid objective functional. Fig. 1b shows the case in which the gradients of the two functionals are equal. In this case, a properly chosen coarse grid objective functional can upper bound the fine grid objective functional, and the coarse grid iteration is ensured to reduce the value of the fine grid objective functional. A V-cycle multigrid algorithm can be easily got by replacing the coarser grid relaxation operator Rg+1 of this two-grid algorithm with another two-grid algorithm recursively. In the inversion process, the V-cycle multigrid method moves from the finest grid to the coarsest grid, and subsequently moves back to the finest grid. The corresponding algorithm is defined as follows: Algorithm 4.1. The V-cycle multigrid algorithm. (1) Initialization. Given q(0) , μ(0) , ud(0 ) = ud , r(0 ) = 0, set g = 0. (g) ). (2) Iterate τ (g) times by the relaxation operator Rg : q(g) ←− Rg (q(g) , K (3) Let q(g+1 ) = Ig

g+1 (g) q ,

(g+1 )

and compute ud

, μ(g+1 ) , r(g+1 ) by Eqs. (4.7), (4.9) and (4.11), respectively.

(g+1 ) ), and go to Step 3 with g = g + 1 (4) Iterate τ (g+1 ) times by the relaxation operator Rg+1 : q(g+1 ) ←− Rg+1 (q(g+1 ) , K until g = G − 1. g+1 (g+1 ) (5) Let q(g+1 ) = Ig q(g) , and compute ud , μ(g+1 ) , r(g+1 ) by Eqs. (4.7), (4.9) and (4.11), respectively. (g+1 ) ). (6) Iterate the relaxation operator Rg+1 until convergence condition is satisfied: q(g+1 ) ←− Rg+1 (q(g+1 ) , K g (7) Carry out the correction: q(g) = q(g) + Ig+1 (q(g+1 ) − q(g+1 ) ). (g) ), and go to Step 7 (8) Iterate the relaxation operator Rg until convergence condition is satisfied: q(g) ←− Rg (q(g) , K with g = g − 1 until g = 0.

T. Liu / Applied Mathematics and Computation 320 (2018) 271–281

277

Fig. 2. The description of the V-cycle multigrid algorithm.

The procedure of Algorithm 4.1 is graphically described by Fig. 2. As shown in Fig. 2, the first half of Algorithm 4.1 (Steps 2–4) is to perform the optimization for a fixed number τ (g) iterations, because the convergence is not required for this stage. Generally speaking, this fixed number τ (g) is usually chosen to be a small value, sometimes even zero. However, the convergence is required for the second half of Algorithm 4.1 (Steps 5–8), that is, the optimization in this stage is performed until the convergence condition is satisfied. 5. Numerical experiments In this section we solve some numerical examples to demonstrate the efficiency and effectiveness of the nonlinear multigrid method. The test problem is Eq. (1.1) with  = [0, 1] × [0, 1], T = 0.1, φ (x ) = sin(π x ) sin(π y ), and s(x, t ) = 0. The nonlinear S-shaped flux functions of Buckley–Leverett type in the x, y directions are chosen as

ϕ (u ) =

u2 (1 − 5(1 − u2 )) u2 + ( 1 − u )2

and

ψ (u ) =

u2

u2 . + ( 1 − u )2

The nonlinear function in the diffusion term is chosen as

N ( u ) = 1 + u ( u − 1 ). The choices of these functions are all from reference [13]. Assume that the initial regularization parameter is μ(0 ) = 10−7 , and take q(0 ) = 5 as the initial guess value. The spatial 1 1 and time step lengths are respectively taken as hx = 48 , hy = 48 , and ht = 0.005. Therefore, there are both 49 grid points in the x direction and in the y direction. In this section, we will consider the two level multigrid method, three level multigrid method, and four level multigrid method, that is G = 1, 2, and 3. The fixed iteration numbers are chosen to be τ ( 0 ) = τ ( 1 ) = τ ( 2 ) = τ ( 3 ) = 1. The observed data ud are generated from the true permeability models by using the finite-difference scheme (2.2). The observed points of data ud are set at the following grid points:

8 8 , 48 ,

48  16 8 , 48 , ..

40 . 8  , , 48 48 48

8 16  , 48 ,

48  16 16 , 48 , ..

40 . 16  , , 48 48 48

... ... .. . ...

8 40  , 48 ,

48  16 40 , 48 , ..

40 . 40  , . 48 48 48

For the selection of these observed points and its relevant application, the interested readers can be referred to [31] and references therein. To illustrate the noise sensitivity, we perform these examples by adding 2% Gaussian noise to the observed data, and then we identify the permeability from the noisy data. The convergence condition is set as

qk+1 − qk  ≤ 10−8 . The relative error is computed by

q(0) − q∗  , q∗ 

278

T. Liu / Applied Mathematics and Computation 320 (2018) 271–281

Fig. 3. The true permeability model in Example 5.1.

Table 1 Comparison of the nonlinear multigrid method and fixed grid method. Example number

Inversion method

Relative error

CPU run time (s)

5.1

Four level multigrid method Three level multigrid method Two level multigrid method Regularized Gauss–Newton method Four level multigrid method Three level multigrid method Two level multigrid method Levenberg Marquardt method

0.0265 0.0265 0.0265 0.0296 0.0641 0.0641 0.0641 0.0752

2995.915 3297.421 4173.523 6979.668 2828.486 3089.876 3802.700 6020.289

5.2

where q(0 ) and q∗ represent, respectively, the inversion result and true permeability model.

Example 5.1. We first consider a level-stratified permeability model containing two interfaces, where the permeability values are as follows:

⎛

5.32 ⎜ .. ⎜ . ⎜5.32 ⎜ ⎜3.98 ⎜ ⎜ .. ⎜ . ⎜ ⎜3.98 ⎜ ⎜4.57 ⎜ . ⎝ .. 4.57

5.32 .. . 5.32 3.98 .. . 3.98 4.57 .. . 4.57

... .. . ... ... .. . ... ... .. . ...

⎞

5.32 ← .. ⎟ . ⎟ ⎟ 5.32⎟← 3.98⎟ ⎟← .. ⎟ . ⎟ ⎟ 3.98⎟← ⎟ 4.57⎟← .. ⎟ ⎠ . 4.57 ←

Line1 .. . Line16 Line17 .. . Line32 Line33 .. . Line48

The graph of this matrix can be plotted by Matlab software and is given in Fig. 3. For comparison, in this example, we investigate the regularized Gauss–Newton method and nonlinear multigrid method (where the relaxation operator is chosen as regularized Gauss–Newton method). Table 1 lists the relative errors and CPU run times of these two methods, and Fig. 4 shows the inversion results by these two methods.

T. Liu / Applied Mathematics and Computation 320 (2018) 271–281

279

Fig. 4. The inversion results in Example 5.1.

Example 5.2. In the second example, we consider a complicated permeability model, where the permeability values are as follows:

⎛

2.77 ⎜ .. ⎜ . ⎜2.77 ⎜ ⎜2.77 ⎜ ⎜ . ⎜ .. ⎜

⎜ ⎜2.77 ⎜ ⎜ .. ⎜ . ⎜ ⎜2.77 ⎜ ⎜2.77 ⎜ . ⎝ .. 2.77 where

⎛

... .. . ... ... .. .

2.77 .. . 2.77 2.77 .. .

2.77 .. . 2.77 6.12 .. .

... .. . ... ... .. . ...

2.77 .. . 2.77 2.77 .. . 2.77

6.12 .. . 6.12 2.77 .. . 2.77

... .. . ... ... .. . .. . .. . ... ... .. . ...

4.09 4.09 4.09 4.09 4.09 4.09 4.09 4.09 4.09 4.09

6.12 6.12 4.09 4.09 4.09 4.09 4.09 4.09 4.09 4.09

6.12 6.12 4.09 4.09 4.09 4.09 4.09 4.09 4.09 4.09

4.09 ⎜4.09 ⎜4.09 ⎜ ⎜4.09 ⎜ ⎜4.09 M=⎜ ⎜4.09 ⎜4.09 ⎜ ⎜4.09 ⎝4.09 4.09

... .. . ... ... .. . M .. . ... ... .. . ...

... .. . ... ... .. . .. . .. . ... ... .. . ...

6.12 6.12 6.12 6.12 4.09 4.09 4.09 4.09 4.09 4.09

2.77 .. . 2.77 6.12 .. .

2.77 .. . 2.77 2.77 .. .

... .. . ... ... .. .

2.77 .. . 2.77 2.77 .. .

6.12 .. . 6.12 2.77 .. . 2.77

2.77 .. . 2.77 2.77 .. . 2.77

... .. . ... ... .. . ...

2.77

6.12 6.12 6.12 6.12 4.09 4.09 4.09 4.09 4.09 4.09

6.12 6.12 6.12 6.12 6.12 6.12 4.09 4.09 4.09 4.09

6.12 6.12 6.12 6.12 6.12 6.12 4.09 4.09 4.09 4.09

⎞

⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ , 2.77 ⎟ ⎟ ⎟ .. ⎟ . ⎟ 2.77 ⎟ ⎟ 2.77 ⎟ ⎟ .. ⎠ .

6.12 6.12 6.12 6.12 6.12 6.12 6.12 6.12 4.09 4.09

⎞

6.12 ← 6.12⎟← 6.12⎟ ⎟← 6.12⎟← ⎟ 6.12⎟← ⎟ 6.12⎟← 6.12⎟ ⎟← 6.12⎟← ⎠ 4.09 ← 4.09 ←

Line21 Line22 Line23 Line24 Line25 Line26 Line27 Line28 Line29 Line30

The graph of this matrix can be plotted by Matlab software and is given in Fig. 5. For comparison, in this example, we investigate the Levenberg Marquardt method and nonlinear multigrid method (where the relaxation operator is chosen as Levenberg Marquardt method). The relative errors and CPU run times of these two methods are also listed in Table 1, and the inversion results by these two methods are shown in Fig. 6. In conclusion, one can state from these two numerical studies that: (i) the nonlinear multigrid method converges faster than the regularized Gauss–Newton method and Levenberg Marquardt method, (ii) the nonlinear multigrid method achieves better solutions than the regularized Gauss–Newton method and Levenberg Marquardt method, (iii) the nonlinear multigrid method is robust with respect to recovering the parameter with the noise in the observation data, and (iv) the nonlinear multigrid method is particularly well suited to the solution of inverse problems because all the common methods (such as regularized Gauss–Newton method and Levenberg Marquardt method) can be employed as the relaxation operator. By

280

T. Liu / Applied Mathematics and Computation 320 (2018) 271–281

Fig. 5. The true permeability model in Example 5.2.

Fig. 6. The inversion results in Example 5.2.

comparing the nonlinear multigrid method with different levels, it can be seen that for the nonlinear multigrid method, the convergence rate increases with the number of levels, and the more levels there are, the more slowly convergence rate increases. Remark 5.1. The spatial grid chosen in this section is dyadic. If the spatial grid is not dyadic, such as cubic, we apply this nonlinear multigrid method by changing the dyadic restriction and prolongation operators to the cubic restriction and prolongation operators. For the cubic spatial grid, the interested readers can be referred to related studies, such as [18]. Remark 5.2. The initial regularization parameter μ(0) is determined by trial and error, i.e., we test several μ(0) values and check which performs best. The relative errors for several values of μ(0) with the nonlinear multigrid method are listed in Table 2. From Table 2, we can see that the initial regularization parameter μ(0 ) = 10−7 performs best.

6. Conclusions This paper presented a nonlinear multigrid method for the inverse problem of the nonlinear convection–diffusion equation in the multiphase porous media flow. Numerical simulations indicate that this nonlinear multigrid method can accelerate convergence, improve the inversion result, and resist noise to a certain degree. It is interesting that the proposed method is very general and can be generalized to more complicated inverse problems.

T. Liu / Applied Mathematics and Computation 320 (2018) 271–281

281

Table 2 Relative errors for several values of μ(0) with the nonlinear multigrid method. Example number 5.1

5.2

μ(0) −7

10 10−6 10−5 10−4 10−3 10−7 10−6 10−5 10−4 10−3

Relative error 0.0265 0.0299 0.0327 0.0411 0.0554 0.0641 0.0775 0.0957 0.1299 0.1627

Acknowledgments This work was supported by the National Natural Science Foundation of China (11271102), the Natural Science Foundation of He-bei Province of China (A2017501001), the Science and Technology Funds for the Universities in He-bei Province of China (Z2015039), the Fundamental Research Funds for the Central Universities (N142303011), the PhD Foundation of Northeast University at Qinhuangdao (XNB2015002). References [1] M. Thullner, M.H. Schroth, J. Zeyer, W. Kinzelbach, Modeling of a microbial growth experiment with bioclogging in a two-dimensional saturated porous media flow field, J. Contam. Hydrol. 70 (2004) 37–62. [2] G.M. Homsy, Viscous fingering in porous media, Annu. Rev. Fluid Mech. 19 (1987) 271–311. [3] V. Kippe, J.E. Aarnes, K.A. Lie, A comparison of multiscale methods for elliptic problems in porous media flow, Computat. Geosci. 12 (2008) 377–398. [4] W.G. Gray, S.M. Hassanizadeh, Macroscale continuum mechanics for multiphase porous-media flow including phases, interfaces common lines and common points, Adv. Water Resour. 21 (1998) 261–281. [5] M.S. Espedal, K.H. Karlsen, Numerical solution of reservoir flow models based on large time step operator splitting algorithm, in: A. Fasano (Ed.), Filtration in Porous Media and Industrial Applications: Lecture Notes in Mathematics, Springer, Berlin, 20 0 0, pp. 9–77. [6] G. Evensen, The ensemble kalman filter: Theoretical formulation and practical implementation, Ocean Dynam. 53 (2003) 343–367. [7] H. Hoel, K.J.H. Law, R. Tempone, Multilevel ensemble kalman filtering, SIAM J. Numer. Anal. 54 (2016) 1813–1839. [8] J. Wang, N. Zabaras, A markov random field model of contamination source identification in porous media flow, Int. J. Heat Mass Transf. 49 (2006) 939–950. [9] R. Kulkarni, A.T. Watson, J.E. Nordtvedt, Estimation of porous media flow functions using NMR imaging data, Magn. Reson. Imaging 16 (1998) 707–709. [10] I. Berre, M. Lien, T. Mannseth, Multi-level parameter structure identification for two-phase porous-media flow problems using flexible representations, Adv. Water Resour. 32 (2009) 1777–1788. [11] H. Krüger, A.A. Grimstad, T. Mannseth, Adaptive multiscale estimation of a spatially dependent diffusion function within porous media flow, Inverse Probl. Sci. Eng. 11 (2003) 341–358. [12] T. Liu, Reconstruction of a permeability field with the wavelet multiscalechomotopy method for a nonlinear convection–diffusion equation, Appl. Math. Comput. 275 (2016) 432–437. [13] T.K. Nilssen, K.H. Karlsen, T. Mannseth, X.C. Tai, Identification of diffusion parameters in a nonlinear convection–diffusion equation using the augmented lagrangian method, Computat. Geosci. 13 (2009) 317–329. [14] A. Jameson, Solution of the euler equations for two dimensional transonic flow by a multigrid method, Appl. Math. Comput. 13 (1983) 327–355. [15] Y. Choi, D. Jeong, J. Kim, A multigrid solution for the cahn-hilliard equation on nonuniform grids, Appl. Math. Comput. 293 (2017) 320–333. [16] M.M. Gupta, J. Zhang, High accuracy multigrid solution of the 3d convection-diffusion equation, Appl. Math. Comput. 113 (20 0 0) 249–274. [17] J.A. Rabi, M.J.S. De Lemos, Optimization of convergence acceleration in multigrid numerical solutions of conductive-convective problems, Appl. Math. Comput. 124 (2001) 215–226. [18] U. Ascher, E. Haber, A multigrid method for distributed parameter estimation problems, ETNA 18 (2003) 1–18. [19] I. Spiliopoulos, D.T. Hristopulos, M.P. Petrakis, A. Chorti, A multigrid method for the estimation of geometric anisotropy in environmental data from sensor networks, Comput. Geosci. 37 (2011) 320–330. [20] A. Mewes, B. Kulessa, J.D.M. Kinley, A.M. Binley, Anisotropic seismic inversion using a multigrid monte carlo approach, Geophys. J. Int. 183 (2010) 267–276. [21] S.F.M. Cormick, J.G. Wade, Multigrid solution of a linearized, regularized least-squares problem in electrical impedance tomography, Inverse Probl. 9 (1993) 697–713. [22] L. Borcea, A nonlinear multigrid for imaging electrical conductivity and permittivity at low frequency, Inverse Probl. 17 (2001) 329–359. [23] J.C. Ye, C.A. Bouman, K.J. Webb, R.P. Millane, Nonlinear multigrid algorithms for bayesian optical diffusion tomography, IEEE T. Image Process. 10 (2001) 909–922. [24] S. Oh, A.B. Milstein, C.A. Bouman, K.J. Webb, A general framework for nonlinear multigrid inversion, IEEE T. Image Process. 14 (2005) 125–140. [25] W. Briggs, V. Henson, S.M. Cormick, A Multigrid Tutorial, second ed., SIAM, Philadelphia, 20 0 0. [26] B. Engquist, S. Osher, One-sided difference approximations for nonlinear conservation laws, Math. Comput. 36 (1981) 321–351. [27] A.N. Tikhonov, V.Y. Arsenin, Solution of Ill-posed Problems, John Wiley and Sons, New York, 1977. [28] H.W. Engl, K. Kunisch, A. Neubauer, Convergence rates for tikhonov regularisation of non-linear ill–posed problems, Inverse Probl. 5 (1989) 523–540. [29] B. Kaltenbacher, Some newton-type methods for the regularization of nonlinear ill–posed problems, Inverse Probl. 13 (1997) 729–753. [30] M. Hanke, A regularizing levenberg-marquardt scheme, with applications to inverse groundwater filtration problems, Inverse Probl. 13 (1997) 79–95. [31] T.K. Nilssen, Parameter estimation with the augmented lagrangian method and a study of some fourth order problems. Doktorarbeit, Universitetet i Bergen 2001.