自伴随Eluer

Comput Visual Sci DOI 10.1007/s00791-007-0078-5 REGULAR ARTICLE A multiscale Eulerian–Lagrangian localized adjoint method for transient advection–diffusion equations with oscillatory coefficients Hong Wang · Yabin Ding · Kaixin Wang · Richard E. Ewing · Yalchin R. Efendiev Received: 31 January 2005 / Accepted: 2 May 2006 © Springer-Verlag 2007 Abstract We develop a multiscale Eulerian–Lagrangian localized adjoint method for transient linear advection– diffusion equations with oscillatory coefficients, which arise in mathematical models for describing flow and transport through heterogeneous porous media, composite material design, and other applications. Keywords Advection–diffusion equations · Eulerian–Lagrangian method · Multiscale method · Reservoir simulation · Porous medium flow 1 Introduction Transient advection–diffusion partial differential equations arise in mathematical models for describing groundwater hydrology, environmental modeling and remediation, petro- leum reservoir simulation and many other fields [2,6,9]. These problems admit solutions with moving steep fronts and complicated structures. Centered finite difference or finite element methods often generate numerical solutions with severe nonphysical oscillations. In industrial applications, Communicated by G. Wittum. H. Wang (B) · Y. Ding Department of Mathematics, University of South Carolina, Columbia, SC 29208, USA e-mail: hwang@math.sc.edu K. Wang School of Mathematics and System Sciences, Shandong University, Jinan, Shandong 250100, China R. E. Ewing · Y. R. Efendiev Institute for Scientific Computation, Texas A&M University, College Station, TX 77843-3404, USA upstream weighting techniques are commonly used to stabilize the numerical approximations to these systems in large-scale simulators. But upwind methods tend to produce excessive numerical diffusion, which smears out moving steep fronts of the solutions to the governing equations and introduces spurious grid orientation effect [6,17]. In addition, advection–diffusion equations for describing porous medium flow through heterogeneous porous media usually contain solutions with multiple spatial and tempo- ral scales. Consequently, a direct numerical solution stra- tegy often requires extremely refined space grids and time steps, due to the wide spectrum of spatial and temporal scales in the solutions and the advection–diffusion feature of the governing equations. The Eulerian–Lagrangian localized adjoint method (ELLAM) [4] provides a general framework for developing characteristic methods to solve transient advection–diffusion equations with general boundary conditions in a mass- conservative manner. The ELLAM schemes have been suc- cessfully applied in the numerical modeling of petroleum reservoir simulation and groundwater contaminant transport and remediation. They generate accurate numerical solutions even if large time steps and spatial grids are used, and are very competitive with many numerical methods [3,8,14–16]. The multiscale finite element (or volume) method (MsFEM) [5,10,11] was introduced in recent years for sol- ving partial differential equations with multiscale solutions. These methods aim at obtaining large scale solutions accu- rately and effectively without resolving small scale details. These ideas were carried out via the construction of basis functions that capture small scale information within each cell, so the effect of small scales on the large scale is captured correctly. We develop a multiscale Eulerian–Lagrangian localized adjoint method (MsELLAM) for a transient linear 123 H. Wang et al. advection–diffusion equation with oscillatory coefficients. The goal of this paper is to expose the fundamental idea in the development of such a scheme within the framework of the ELLAM formulation and the MsFEM, so we focus on a simple model problem in most of the paper. At the end of this paper we briefly discuss extensions of the MsELLAM scheme for coupled systems modeling porous medium flow. 2 An ELLAM formulation In this section we revisit the idea of the ELLAM framework [4] in the context of one-dimensional advection–diffusion equations with oscillatory coefficients ct + ( V ( x, x ε ) c − D ( x, x ε ) cx ) x = f, x ∈ (a, b), t ∈ (0, T ], (1) where V is velocity field, D is the diffusion coefficient, f is a given source or sink term, and c(x, t) is a measure of concen- tration of a dissolved substance. Due to the heterogeneity of the porous medium, the velocity V and diffusion coefficient D may be random or highly oscillatory. Consequently, the solution could contain multiple scales. In the model problem (1), we assume the problem to have two scales: a scale of O(1) that represents the normal scale, and a scale of O(ε) that represents a fast changing scale. Note that we do not impose a periodicity assumption on the problem. The ELLAM framework can treat advection–diffusion equations with general inflow and outflow boundary condi- tions in a mass-conservative manner [4,14]. Since the goal of this paper is to explore how to utilize the idea of the multiscale finite element method [10] within the ELLAM framework to design a multiscale ELLAM scheme, we skip the treatment of boundary conditions by, e.g., assuming an initial condition with compact support or periodic boundary conditions. Equation (1) is also subject to the initial condition c(x, 0) = c0(x), x ∈ [a, b]. (2) 2.1 Weak formulation We begin by defining a space–time partition: xi := a + i∆x, 0 ≤ i ≤ I, ∆x := b − aI tn := n∆t, 0 ≤ n ≤ N , ∆t := TN . (3) Due to various mathematical, numerical, and computational constraints, the size of spatial grids and time steps ∆x and ∆t often has to be chosen much larger than the fine scale ε in practical simulations. We multiply Eq. (1) by space–time test functions w that vanish outside [a, b] × (tn−1, tn] and are discontinuous in time at time step tn−1 and integrate the resulting equation to get tn∫ tn−1 b∫ a ( ct + ( V ( x, x ε ) c − D ( x, x ε ) cx ) x ) wdxdt = tn∫ tn−1 b∫ a f wdxdt. (4) Integration by parts in both space and time yields a weak formulation b∫ a c(x, tn)w(x, tn)dx + tn∫ tn−1 b∫ a D ( x, x ε ) cxwx dxdt − tn∫ tn−1 b∫ a c ( wt + V ( x, x ε ) wx ) dxdt = b∫ a c(x, tn−1)w(x, t+n−1)dx + tn∫ tn−1 b∫ a f wdxdt. (5) Here w(x, t+n−1) := limt→t+n−1 w(x, t) accounts for the dis- continuity of w(x, t) at time tn−1. Based on the idea of the localized adjoint method or opti- mal test function method [4,12], the test functions w should be chosen from the solution space of the homogeneous adjoint equation of Eq. (1) − wt − V ( x, x ε ) wx − ( D ( x, x ε ) wx ) x = 0. (6) 2.2 Operator splitting Because the solution space of ordinary differential equa- tions is finite-dimensional, the development of optimal test function methods for the boundary-value problem of one- dimensional ordinary differential equations is relatively simple [4,12]. In contrast, the solution space of the partial dif- ferential equation (6) is infinite-dimensional. Since the objec- tive of a numerical procedure is to derive a finite-dimensional approximation, only a finite number of test functions should be chosen. Different choices of test functions lead to different classes of approximations. 2.2.1 An optimal test function operator-splitting In the context of the adjoint equation (6), one option is to split the adjoint operator into a spatial operator and a tempo- ral operator. Namely, we require the test functions w(x, t) to satisfy the following system that consists of two sub- equations within each element (xi−1, xi ) × (tn−1, tn] 123 A multiscale ELLAM scheme (i = 1, 2, . . . , I ) wt = 0, −V ( x, x ε ) wx − ( D ( x, x ε ) wx ) x = 0. (7) The second equation in (7) is a singularly-perturbed, two- scale, steady-state advection–diffusion equation. It defines the spatial configuration of the test functions w. This type of splitting leads to a class of optimal test function methods involving exponential upstream weighting in space. The cor- responding test functions w typically exhibit exponential boundary layers near the inter-element boundary, which require a fine grid for reasonable resolution. The first equation in (7) defines the temporal variation of the test functions w. Numerical methods resulting from these test functions tend to have large temporal truncation errors, numerical dispersion, and phase errors [4,14,15]. The combination of multiscale feature and advection dominance of the problems introduces extra numerical difficulties. 2.2.2 An Eulerian–Lagrangian operator-splitting An alternative operator splitting of the adjoint equation (6) is to split the adjoint operator into a first-order hyperbolic operator and a second-order diffusion operator. Namely, we require the test functions w(x, t) to satisfy the following sys- tem that consists of two sub-equations within each Eulerian– Lagrangian element Ωi , which is a space–time strip emana- ting backward from (xi−1, xi ) along the characteristics from time step tn to time step tn−1 −wt − V ( x, x ε ) wx = 0, − ( D ( x, x ε ) wx ) x = 0. (8) In the ELLAM framework, the equations in (8) are impo- sed locally within the interior of each element Ωi . As in the case of optimal test function splitting, the second equation in (8) defines the spatial configuration of the test functions w. This equation is now an elliptic steady-state diffusion equa- tion, in contrast to a steady-state advection–diffusion equa- tion. Consequently, the test functions w are not expected to exhibit boundary layers near the inter-element boundary. The first equation in (8) defines the temporal variation of the test functions w, which implies the Lagrangian feature of the test functions. We note that in the case of scale separation (e.g., problems with periodic coefficients) the solution of the local problems (7) and (8) can be approximated using auxi- liary problems, which arise in homogenization. This approxi- mation will reduce the computational cost of the proposed method. 3 An MsELLAM tracking fine-scale velocity The objective of this section is to derive an MsELLAM for problem (1), (2) by incorporating multiscale-handling capa- bility into the ELLAM framework. 3.1 Definition of test functions One of the key issues in the development of MsELLAM is the choice of test functions in the ELLAM framework. In order for the MsELLAM to handle multiscale advection– diffusion equations, we utilize the idea of multiscale finite element method (MsFEM) [10] in the construction of the test functions in the MsELLAM. 3.1.1 Spatial configuration of the test functions The second equation in the Eulerian–Lagrangian operator- splitting (8), which defines spatial configuration of the test functions w in the ELLAM framework, is a self-adjoint ellip- tic equation. The original MsFEM was developed for such an equation [10]. We naturally follow the idea of MsFEM to define the spatial configuration of the test functions w at time step tn . Because the mesh size of the coarse grid ∆x > ε, the test functions in a conventional finite element method, which is linear on [xi−1, xi ], cannot handle fine-scale behavior of the solutions. In the MsFEM a subdivision xi, j := xi−1 + j∆x f , 0 ≤ j ≤ J, ∆x f := ∆xJ (9) is introduced in [xi−1, xi ], where ∆x f = O(ε) is chosen to account for the fine-scale behavior. Then two piecewise- linear basis functions w(1)i (x, tn) and w (2) i (x, tn) are deter- mined by using a finite element method to solve the second equation in (8) in [xi−1, xi ] with respect to the subdivision (9) and the boundary condition w (1) i (xi−1, tn) = 1, w(1)i (xi , tn) = 0, w (2) i (xi−1, tn) = 0, w(2)i (xi , tn) = 1. (10) The shape functions w(1)i (x, tn) and w (2) i (x, tn) are restricted to the element [xi−1, xi ], and are extended outside of the element by 0. These shape functions for i = 0, 1, . . . , I form a basis for the MsFEM with the degrees of freedom being still located at coarse spatial grid nodes xi for i = 0, 1, . . . , I . In this way, the basis functions naturally build the multiscale behavior into their construction and reflect the multiscale behavior due to diffusion. 3.1.2 Temporal variation of the test functions Once the spatial configuration of the test functions w(x, tn) is determined at time step tn , the first equation in the Eulerian– 123 H. Wang et al. Lagrangian splitting (8) defines the temporal variation of w(x, t) from time step tn to time step tn−1 by extending them backward along the characteristic curves y = r(t; x, tn) defi- ned by dy dt = V ( y, y ε ) , y ∣∣∣ t=tn = x . (11) In the case of multiscale transient advection–diffusion equations, the velocity field exhibits multiscale or oscillatory behavior. Consequently, the evaluation of the test functions w, which is carried out via a characteristic tracking, requires an accurate resolution of problem (11) determined by the fine-scale velocity field V . Since the size of the global time step ∆t > ε, we introduce a local time-step partition tn,k := tn−1 + k∆t f , 0 ≤ k ≤ K , ∆t f := ∆tK (12) where the size of the local time step ∆t f = O(ε) is chosen to handle the multiscale behavior in time. 3.2 Numerical scheme With the test functions defined in Sect. 3.1, we now derive an MsELLAM scheme. For clarity of exposition, we use the variable y to represent the spatial coordinate of any point in the domain. We use the change of variable y = r(θ; x, tn) to evaluate the source term in Eq. (5) to obtain ∫ ∫ Ωi f (y, t)w(y, t)dydt = tn∫ tn−1 r(t;xi ,tn)∫ r(t;xi−1,tn) f (y, t)w(y, t)dydt = K∑ k=1 tn,k∫ tn,k−1 r(t;xi ,tn)∫ r(t;xi−1,tn) f (y, t)w(y, t)dydt. (13) Then we can use a numerical quadrature to discretize the temporal integral on each local time interval [tn,k−1, tn,k] for k = 1, 2, . . . , K . We handle the diffusion term in a similar manner ∫ ∫ Ωi D ( y, y ε ) cy(y, t)wy(y, t)dydt = tn∫ tn−1 r(t;xi ,tn)∫ r(t;xi−1,tn) D ( y, y ε ) cy(y, t)wy(y, t)dydt = K∑ k=1 tn,k∫ tn,k−1 r(t;xi ,tn)∫ r(t;xi−1,tn) D ( y, y ε ) cy(y, t)wy(y, t)dydt. (14) Because the test functions w satisfy the first equation of (8) (at least approximately), the last term on the left-hand side of (5) vanishes (or at least is the same order as the global truncation error of the numerical scheme [7,13]. Hence, we drop this term in the numerical scheme. We then obtain an MsELLAM scheme by incorporating Eqs. (13) and (14) into Eq. (5). 4 An MsELLAM tracking coarse-scale velocity The MsELLAM scheme developed in the previous section requires a characteristic tracking of the fine-scale velocity field. In this section we explore the possibility of developing an MsELLAM for problem (1), (2), which requires only tra- cking the coarse-grid velocity field. 4.1 Definition of test functions One of the key issues in the development of MsELLAM is the choice of test functions in the ELLAM framework, which in turn is the result of an appropriate operator-splitting. 4.1.1 An alternative Eulerian–Lagrangian operator-splitting To develop an MsELLAM scheme that requires only tracking the coarse-grid velocity field, we decompose the fine-scale velocity field V as V = V¯ + V˜ . (15) Here V¯ and V˜ represent the global, smooth part and the local, oscillatory part of the velocity field V , respectively. We revise the Eulerian–Lagrangian operator-splitting (8) as follows −wt − V¯ wx = 0, −V˜x − (Dwx )x = 0. (16) In the ELLAM framework, the equations in (8) are imposed locally within the interior of each element Ωi . We note that a somewhat related modified Eulerian–Lagrangian operator- splitting was previously used in [7] for the purpose of improving the accuracy of the numerical approximation and simplifying the characteristic tracking algorithm. 4.1.2 Spatial configuration of the test functions We note that the second equation in the operator-splitting (16) is now a steady-state advection–diffusion equation. We use the MsFEM idea to define the spatial configuration of the basis functions w at time step tn . In the current context, we solve the second equation in (16) by a standard finite ele- ment method with the fine-scale subdivision (9) on interval 123 A multiscale ELLAM scheme [xi−1, xi ]. Because the fine-scale grid size ∆x f = O(ε) is chosen to be small enough so that the grid Peclet number is of O(1), we do not expect that the basis functions exhibit nonphysical oscillations [12]. 4.1.3 Temporal variation of the test functions Once the spatial configuration of the test functions w(x, tn) is determined at time step tn , the first equation in the modi- fied Eulerian–Lagrangian splitting (16) defines the tempo- ral variation of w(x, t) from time step tn to time step tn−1 by extending them backward along the characteristic curves y = r(t; x, tn) defined by dy dt = V¯ ( y, y ε ) , y ∣∣∣ t=tn = x . (17) 4.2 Numerical scheme The diffusion term and source term in the reference equa- tion (5) can be handled in a similar manner to Eqs. (13) and (14). The first equation in the modified Eulerian–Lagrangian operator-splitting (16) implies tn∫ tn−1 b∫ a c(wt + V¯ wx )dxdt (18) vanishes, or at least, is within the global truncation error of the numerical scheme [7]. The oscillatory part V˜ of the velocity field on the left-hand side of the reference equation (5) remains, which can be approximated as follows ∫ ∫ Ωi V˜ ( y, y ε ) c(y, t)wy(y, t)dydt = tn∫ tn−1 r(t;xi ,tn)∫ r(t;xi−1,tn) V˜ ( y, y ε ) c(y, t)wy(y, t)dydt = K∑ k=1 tn,k∫ tn,k−1 r(t;xi ,tn)∫ r(t;xi−1,tn) V˜ ( y, y ε ) c(y, t)wy(y, t)dydt. (19) We handle the temporal integral in this term as we did with the source term in Eq. (13) and with the diffusion term in (14) in Sect. 3. 5 Numerical experiments We simulate the transport of a one-dimensional Gaussian pulse over the spatial domain [0, 4]. In the example runs the time interval is [0, T ] = [0, 1]. The initial condition is given by c0(x) = exp ( − (x − x0) 2 2σ 2 ) , (20) where the center x0 = 0.5 and the spread σ = 0.1 that determines the steepness of the Gaussian pulse. A two-scale oscillatory velocity filed V (x/ε) and diffu- sion coefficient D(x/ε) is chosen as follows V ( x ε ) = 1 + α cos (2πx ε ) , D ( x ε ) = 0.01 1 + 0.8 sin (2πx ε ) , (21) where α = 0.2 which determines the magnitude of the oscil- latory velocity field. We use ε = 0.02, 0.01 and 0.005 to investigate the convergence rate. The source term f is cho- sen to be zero. 5.1 Results with ELLAM We use an ELLAM scheme with extremely refined spatial grid and time step to generate a reference solution, which is presented in Fig. 1. In Table 1 and Fig. 1 we present numeri- cal results generated by an ELLAM scheme with an accurate characteristics tracking using a second-order Runge–Kutta formula. These numerical results show that the ELLAM scheme converges only when the meshsize ∆x < ε. This is fully understandable and can be explained as follows. Although ELLAM schemes have been successfully applied to simulate transient advection–diffusion equations in dif- ferent applications [3,8,14–16], they are not designed to solve advection–diffusion equations with multiple scales. This is similar to the case why conventional finite element methods do not work well for multiscale elliptic equations. 0.8 1 1.2 1.4 1.6 1.8 2 2.2 0 0.1 0.2 0.3 0.4 0.5 Reference Solution dx=1/5 dx=1/10 dx=2/25 dx=1/20 dx=1/25 Fig. 1 Results by an ELLAM scheme with different grids 123 H. Wang et al. Table 1 Results with an ELLAM (ε = 1100 , ∆t = 11000 ) ∆x L1 error L2 error L∞ error 1 5 0.0809 0.0784 0.1758 4 25 0.0888 0.0908 0.1981 1 10 0.0651 0.0753 0.1432 2 25 0.0518 0.0602 0.1202 1 20 0.0499 0.0564 0.1225 1 25 0.1269 0.1390 0.2263 1 50 0.0346 0.0390 0.0729 1 100 0.0581 0.0644 0.1376 2 250 0.0484 0.0551 0.1229 1 20