Multirate partially explicit scheme for multiscale flow problems

For time-dependent problems with high-contrast multiscale coefficients, the time step size for explicit methods is affected by the magnitude of the coefficient parameter. With a suitable construction of multiscale space, one can achieve a stable temporal splitting scheme where the time step size is independent of the contrast. Consider the parabolic equation with heterogeneous diffusion parameter, the flow rates vary significantly in different regions due to the high-contrast features of the diffusivity. In this work, we aim to introduce a multirate partially explicit splitting scheme to achieve efficient simulation with the desired accuracy. We first design multiscale subspaces to handle flow with different speed. For the fast flow, we obtain a low-dimensional subspace with respect to the high-diffusive component and adopt an implicit time discretization scheme. The other multiscale subspace will take care of the slow flow, and the corresponding degrees of freedom are treated explicitly. Then a multirate time stepping is introduced for the two parts. The stability of the multirate methods is analyzed for the partially explicit scheme. Moreover, we derive local error estimators corresponding to the two components of the solutions and provide an upper bound of the errors. An adaptive local temporal refinement framework is then proposed to achieve higher computational efficiency. Several numerical tests are presented to demonstrate the performance of the proposed method.

Authors

• 10 publications
• 8 publications
• Contrast-independent partially explicit time discretizations for multiscale wave problems

In this work, we design and investigate contrast-independent partially e...
02/25/2021 ∙ by Eric T. Chung, et al. ∙ 0

• Partially Explicit Time Discretization for Time Fractional Diffusion Equation

Time fractional PDEs have been used in many applications for modeling an...
08/30/2021 ∙ by Jiuhua Hu, et al. ∙ 0

• Temporal Splitting algorithms for non-stationary multiscale problems

In this paper, we study temporal splitting algorithms for multiscale pro...
11/11/2020 ∙ by Yalchin Efendiev, et al. ∙ 0

• Contrast-independent partially explicit time discretizations for nonlinear multiscale problems

This work continues a line of works on developing partially explicit met...
08/28/2021 ∙ by Eric T. Chung, et al. ∙ 0

• Design of High-Order Decoupled Multirate GARK Schemes

Multirate time integration methods apply different step sizes to resolve...
04/20/2018 ∙ by Arash Sarshar, et al. ∙ 0

• Upscaling errors in Heterogeneous Multiscale Methods for the Landau-Lifshitz equation

In this paper, we consider several possible ways to set up Heterogeneous...
04/07/2021 ∙ by Lena Leitenmaier, et al. ∙ 0

• Wavelet-based Edge Multiscale Parareal Algorithm for Parabolic Equations with Heterogeneous Coefficients

We propose in this paper the Wavelet-based Edge Multiscale Parareal (WEM...
03/23/2020 ∙ by Guanglian Li, et al. ∙ 0

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

Modeling of flow and transport in complicated porous media in various physical and engineering applications encounters problems with multiscale features. In particular, the properties of the underlying media, such as thermal diffusivity or hydraulic conductivity, have values across different magnitudes. This poses challenges in the numerical simulation since the high contrast feature of the heterogeneous media introduces stiffness for the system. In terms of temporal discretization, the time-stepping depending on the magnitude of the multiscale coefficient is needed for explicit schemes. For the spatial discretization, many multiscale methods [27, 25, 10, 22, 5, 20] are introduced to handle the issue. The multiscale model reduction methods include both local [22, 2, 1, 4, 21] and global [26, 7, 9, 8] approaches to reduce computational expenses. The idea is to construct reduced order models to approximate the full fine-scale model and achieve efficient computation. Among these methodologies, the family of generalized multiscale finite element methods (GMsFEM) [20, 13, 16, 17]

are proposed to effectively address multiscale problem with high-contrast parameters. It first formulates some local problems on coarse grid regions to get snapshot basis that can capture the heterogeneous properties, and then designs appropriate spectral problems to get important modes in the snapshot space. The GMsFEM approach share some similarities with multi-continuum methods. The basis functions can recognize the high-contrast features such as channels that need to be represented individually. The convergence of the GMsFEM depends on the eigenvalue decay, and the small eigenvalues correspond to the high permeable channels.

To construct multiscale basis such that the convergence of the method is independent of the contrast and linearly decreases with respect to mesh size under suitable assumptions, the constraint energy minimizing GMsFEM (CEM-GMsFEM) was initiated[14, 11]. This approach begins with a suitable choice of auxiliary space, where some local spectral problems in coarse blocks are solved. The auxiliary space includes the minimal number of basis functions to identify the essential information of the channelized media. Then it will be used to compute the solutions of constraint energy minimizing problem in some oversampling coarse regions to handle the non-decaying property. The resulting localized solutions form the multiscale space.

To adapt the CEM-GMsFEM for flow-based upscaling, the nonlocal multicontinuum upscaling method (NLMC) [12] is proposed by modifying the above framework. The idea is to use simplified auxiliary space by assuming that each separate fracture network within a coarse grid block is known. The auxiliary basis are piecewise constants corresponding to fracture networks and matrix, which are called continua. Then the local problems are formulated for each continuum by minimizing the local energy subject to appropriate constraints. This construction returns localized basis functions which can automatically identify each continuum. Further, due to the property of the NLMC basis, this approach will provide non-local transmissibilities which describe the transfer among coarse blocks in an oversampled region and among different continua.

Consider the time-dependent problem with high-contrast coefficients, there have been various approaches to handle multiscale stiff systems [3, 6, 23, 24]. Recently, a temporal splitting method is combined with the spatial multiscale method [15] to produce a contrast-independent partially explicit time discretization scheme. It splits the solution of the problem into two subspaces which can be computed using implicit and explicit methods, instead of splitting the operator of the equation directly based on physics [32, 33, 34, 28, 29]. The multiscale subspaces are carefully constructed. The dominant basis functions stem from CEM-GMsFEM which have very few degrees of freedom and are treated implicitly. The additional space as a complement will be treated explicitly. It was shown that with the designed spaces, the proposed implicit-explicit scheme is unconditionally stable in the sense that the time step size is independent of the contrast. Following a similar idea in [15], in this work, we will propose a multirate time-stepping method for the multiscale flow problem.

Multirate time integration method has been studied extensively for the flow problems [18, 19, 30, 31]. Based on different splittings of the target equation, multiple time stepping are utilized in different parts of the spatial domain according to computational cost or complexity of the physics.

In this paper, we integrate the multirate approach with multiscale space construction. Due to the high contrast property of the coefficients, the solutions pass through different regions of the porous medium with different speeds in the flow problem. Different from the previous approach [15], where the multiscale basis functions are formulated for dominant features and missing information, we propose to design multiscale spaces to handle the fast and slow components of the slow separately. We use the simplified auxiliary space containing piecewise constant functions as in the NLMC framework. We only keep the basis representing the high-diffusive component in the first space and adopt an implicit time discretization scheme. The second space consists of bases representing the remaining component, it will take care of the slow flow and the corresponding degrees of freedom are solved explicitly. Next, we introduce a multirate approach where different time step sizes are employed in the partially explicit splitting scheme, such that different parts of the solution are sought with time steps in line with the dynamics. We start with a coarse step size for both equations and refine local coarse time blocks based on some error estimators. With a finer discretization, the accuracy of the approximation can be improved. We analyze the stability of the multirate methods for all four cases when we use coarse or fine time step size alternatively for the implicit and explicit parts of the splitting scheme. It shows that the scheme is stable as long as the coarse time step size satisfies some suitable conditions independent of the contrast. Moreover, we propose an adaptive algorithm for the implicit-explicit scheme by deriving error estimators based on the residuals. The two error estimators corresponding to the two components of the solutions can provide an upper bound of the errors. Compared with uniform refinement, an adaptive refining algorithm can enhance the efficiency significantly. Several numerical examples are presented to demonstrate the effectiveness of the proposed adaptive method.

The paper is organized as follows. In Section 2, we describe the problem setup and the partially explicit scheme. The construction of the multiscale spaces is discussed in Section 3. In Section 4, the multirate method is presented, the subsection 4.1 is devoted to the stability analysis and the subsection 4.3 presents the adaptive algorithm. Numerical tests are shown in Section 5. A conclusion is drawn in Section 6.

2 Problem Setup

Consider the parabolic equation

 dudt−∇⋅(κ∇u) =f on Ω×(0,T] u =0 on ∂Ω×(0,T] u =u0 on ∂Ω×{0}

where is a heterogeneous coefficient with high contrast, that is, the value of the conductivity/permeability in different regions of can differ in magnitudes.

For a finite element space , the discretization in space leads to seeking such that

 ddt(u,v)+a(u,v) =(f,v),∀v∈V,t∈(0,T] u(0,⋅) =u0

where .

Now consider a coarse spatial partition of the computational domain , we will construct suitable multiscale basis functions on and form a multiscale space which is an approximation of the finite element space. Let be the time step size. The discretization in the space with implicit backward Euler scheme in time reads

 (uk+1H−ukHτ,v)+a(uk+1H,v)=(fk+1,v),∀v∈VH (1)

where is the number of time steps, and . It is well-known that this implicit scheme is unconditionally stable.

Suppose the multiscale space can be decomposed into two subspaces

 VH=VH,1+VH,2,

then a partial explicit temporal splitting scheme [15] is to find and , for all satisfying

 ⎛⎝uk+1H,1−ukH,1τ,v1⎞⎠+⎛⎝ukH,2−uk−1H,2τ,v1⎞⎠+a(uk+1H,1+ukH,2,v1)=(fk+1,v1), (2)
 ⎛⎝uk+1H,2−ukH,2τ,v2⎞⎠+⎛⎝ukH,1−uk−1H,1τ,v2⎞⎠+a((1−ω)ukH,1+ ωuk+1H,1+ukH,2,v2) (3) =(fk+1,v2),

, where is a customized parameter. In the case , the two equations are decoupled, and can be solved simultaneously. In the case , the second equation depends on the solution , thus the two equations will be solved sequentially.

The solution at time step will be . It was shown in [15] that under appropriate choices of the multiscale spaces and , the above implicit-explicit scheme resulted from the temporal splitting method for multiscale problems are stable with time step independent of contrast. In [15], the dimension of is low and it contains some dominant multiscale basis functions, the second space includes additional bases representing the missing information. In this paper, we will construct multiscale spaces corresponding to different time scales, where the fast and slow parts of the solution are treated separately.

3 Construction of multiscale spaces

In this section, we will present the construction of multiscale spaces.

3.1 Multiscale basis for Vh,1

We will first discuss the basis construction for based on the contraint energy minimizing GMsFEM (CEM-GMsFEM) [14] and the nonlocal multicontinuum method (NLMC)[12].

To start with, we introduce some notations for the fine and coarse discretization of the computational domain . Let be a coarse partition with mesh size , and be a conforming refinement of with mesh size , where . Denote by () the set of coarse blocks in , and is an oversampled region with respect to each , where the oversampling part contains a few layers of coarse blocks neighboring . Let be the restriction of on .

Under the framework of CEM-GMsFEM, one first constructs an auxiliary space. Consider the spectral problem

 ai(ϕ(i)aux,k,v)=λiksi(ϕ(i)aux,k,v),∀v∈V(Ki), (4)

where and are corresponding eigenpairs, and

 ai(u,v)=∫Ki∇u⋅∇v,si(u,v)=∫Ki~κuv,

with , and denotes the multiscale partition of unity function. Upon solving the spectral problem, we arrange the eigenvalues of (4) in an ascending order, and select the first eigenfunctions to form the auxiliary basis functions. Define , where and is the number of coarse elements. The the global auxiliary space . We note that the auxiliary space needs to be chosen appropriately in order to get good approximation results. That is, the first few basis functions corresponding to small eigenvalues (representing all the channels) have to be included in the space.

Define a projection operator satisfying

 πi(u)=li∑k=1si(u,ϕ(i)aux,k)si(ϕ(i)aux,k,ϕ(i)aux,k)ϕ(i)aux,k,∀u∈V.

Then we let , and . Define the null space of to be :

 ~V={v∈V|π(v)=0}.

Let the global basis be the solution of the optimization problem

 ψ(i)glo,j=argmin{a(v,v) |v∈V0(K+i),s(v,ϕ(i)aux,k)=1 and s(v,ϕ(i′)aux,k′)=0∀i′≠i,k′≠k}.

Define . It can be seen that is -orthogonal to , that is

 a(ψ(i)glo,j,v)=0;∀v∈~V.

Then the CEM multiscale basis is a localization of , and are also computed using the auxiliary space . The idea is to solve the constraint energy minimization problem in a localized region

 a(ψ(i)cem,j,w)+s(w,μ(i)j) =0,∀w∈V(K+i), (5) s(ψ(i)cem,j,ν) =s(ϕ(i)aux,j,ν),∀ν∈V(ℓ)aux% ,

where is an auxiliary basis.

The multiscale space is then , it is an approximation to the global space .

Note that the construction of CEM basis which we have presented here is general and can handle complex heterogeneous permeability field with high contrast. In a fractured media where the media configurations are explicitly known, here the assumption is valid in many real applications, we can consider the simplified construction.

The domain for the media with fracture networks can be represented as follows

 Ω=Ωms⨁l=1dlΩf,l

where the subscripts and denote the matrix and fractures correspondingly. In the fracture regions , the scalar denotes the aperture, and is the number of discrete fracture networks. The permeabilities of matrix and fractures usually have high contrast.

In this setting, some simplified auxiliary basis can be adopted, the constraint energy minimizing basis can be constructed via NLMC [12] and the resulting basis functions which can separate the continua such as matrix and fracture automatically. To be specific, for a given coarse block, we use constants for each separate fracture network, and then a constant for the matrix for the simplified auxiliary space.

Consider an oversampled region , for a coarse element , let be the set of discrete fractures, and be the number of elements in . The NLMC basis are obtained by solving the following localized constraint energy minimizing problem on the fine grid

 a(ψ(i)m,v)+∑Kj⊂K+i⎛⎝μ(j)0∫Kjv+∑1≤n≤mjμ(j)n∫f(j)nv⎞⎠=0, ∀v∈V0(K+i), (6) ∫Kjψ(i)m=δijδm0, ∀Kj⊂K+i, ∫f(j)nψ(i)m=δijδmn, ∀f(j)n∈Fj,∀Kj⊂K+i.

We remark that the resulting basis separates the matrix and fractures automatically, and have spatial decay property. The NLMC basis for fractured media is then .

We denote the the average of the NLMC basis representing the fractures by

 ¯ψ:=1LNc∑i=1mi∑m=1ψ(i)m (7)

, where . Note that in this double summation, the subscript starts from , which indicates the basis functions corresponding to the matrix are not included here. Let . To simplify the notation, we omit the double scripts in and denote the set of basis by

Finally, we define the space as follows:

 VH,1=span{~ψk,1≤k≤L−1}, (8)

where we take away the last basis to remove linear dependency. We remark that, contains basis representing the high contrast regions only. The basis functions corresponding to the matrix and the basis will be included in the second subspace in the following section.

3.2 Multiscale basis for Vh,2

In this section, we will present basis construction for . We first include which are computed from equations (6), and the average basis obtained from equation (7) in . Additionally, we can also conduct some spectral decomposition in the local coarse region to construct more basis corresponding to the information in the non-high-contrast regions. Here be a coarse neighborhood with respect to the -th coarse node. The basis in the spaces needs to satisfy the stability condition (16) and (20). The additional basis are obtained through finding the eigenvalues and corresponding eigenfunctions from the spectral problem

 ∫ωiκ∇z(i)k⋅∇v=η(i)kH2∫ωiz(i)kv (9)

for all .

Again, rearranging eigenvalues in an ascending order, and select the first eigenfunctions correspondingly, we get

 {z(i)k,0≤k≤Li,1≤i≤N}.

The space is then defined to be

 VH,2=span{ψ(i)0,¯ψ,z(i)k,0≤k≤Li,1≤i≤N}. (10)

We remark that, for simplicity, we only use for the second space in our numerical examples.

4 Multirate time stepping for partially explicit scheme

Based on the multiscale spaces constructed in Section 3, we introduce a multirate time stepping partially explicit temporal splitting scheme. Consider the coarse time step size and fine time step size , where . Denote by the fine partition of the time domain by

 0=t0

The coarse partition of the time domain is formed by

 0=T0

Further, we write each coarse time interval where .

The multirate scheme is then defined as follows. For each coarse interval , we are seeking for with a given , such that defined in one of the following four cases which are using coarse time step size in both (2) and (3) (coarse-coarse), using coarse time step size in (2) and using fine time step size in (3) (coarse-fine), using coarse time step size in (2) and using fine time step size in (3) (fine-coarse), using fine time step size in (2) and using fine time step size in both (2) and (3) (fine-fine).

• Case 1 (coarse-coarse): Coarse time step size for (2), coarse time step size for (3). That is, take in both equations, let :

 ⎛⎝unk+1H,1−unkH,1ΔT,v1⎞⎠+⎛⎝unkH,2−unk−1H,2ΔT,v1⎞⎠+a(unk+1H,1+unkH,2,v1)=0, (11) ⎛⎝unk+1H,2−unkH,2ΔT,v2⎞⎠+⎛⎝unkH,1−unk−1H,1ΔT,v2⎞⎠+a(¯unk+1H,1+unH,2,v2)=0,

.

• Case 2 (coarse-fine): Coarse time step size for (2), fine time step size for (3), let :

 ⎛⎝unk+1H,1−unkH,1ΔT,v1⎞⎠+⎛⎝unkH,2−unk−1H,2ΔT,v1⎞⎠+a(unk+1H,1+unkH,2,v1)=0, (12) ⎛⎝un+1H,2−unH,2Δt,v2⎞⎠+⎛⎝unkH,1−unk−1H,1ΔT,v2⎞⎠+a(¯unk+1H,1+unH,2,v2)=0,

, and for .

• Case 3 (fine-coarse): Fine time step size for (2), coarse time step size for (3)

 ⎛⎝un+1H,1−unH,1Δt,v1⎞⎠+⎛⎝unkH,2−unk−1H,2ΔT,v1⎞⎠+a(un+1H,1+unkH,2,v1)=0, (13) ⎛⎝unk+1H,2−unkH,2ΔT,v2⎞⎠+⎛⎝unkH,1−unk−1H,1ΔT,v2⎞⎠+a(¯unk+1H,1+unkH,2,v2)=0,

, and for .

• Case 4 (fine-fine): Fine time step for (2), fine time step for (3). That is, take in both equations, let :

 ⎛⎝un+1H,1−unH,1Δt,v1⎞⎠+⎛⎝unH,2−un−1H,2ΔT,v1⎞⎠+a(un+1H,1+unH,2,v1)=0, (14) ⎛⎝un+1H,2−unH,2ΔT,v2⎞⎠+⎛⎝unH,1−un−1H,1ΔT,v2⎞⎠+a(¯un+1H,1+unH,2,v2)=0,

, and for

We remark that

may not be defined in some of the fine time steps. In this case, we use the linear interpolation of the nearest two coarse time step solutions

to define intermediate time step solutions for .

4.1 Stability for different cases

Consider a coarse time block , the stability of the multirate method for the above mentioned four cases is proved in this subsection.

Let be a constant such that

 γ=supv1∈VH,1,v2∈VH,2(v1,v2)∥v1∥∥v2∥<1 (15)

For case 1 and case 4 defined in section 4.3, following a similar proof in [15], the partially explicit scheme (2)-(3) is stable if

 τsupv∈VH,2∥v∥2a∥v∥2≤1−γ22−ω,

and for case 1, for case 4.

We will show the stability for case 2 and 3 in the following.

4.1.1 Stability for case 2

Use the coarse time step for and use the fine time step for ,

Lemma 1.

The multirate partially explicit scheme in (12) satisfies the stability estimate

 γ2ΔT22∑j=1∥unk+1H,j−unkH,jΔT∥2+12∥unk+1H∥2a≤γ2ΔT22∑j=1∥unkH,j−unk−1H,jΔT∥2+12∥unkH∥2a.

if

 ΔTsupv∈VH,2∥v∥2a∥v∥2≤(1−γ2)mm+1−mω. (16)
Proof.

The equations in (12) can be written as

 (unk+1H,1−unkH,1+unkH,2−unk−1H,2,v1)=−ΔTa(unk+1H,1+unkH,2,v1), (17)
 (m(un+1H,2−unH,2)+unkH,1−unk−1H,1,v2)=−ΔTa((1−ω)unkH,1+ωunk+1H,1+unH,2,v2). (18)

Take in (17), take in (18) and sum over . Then for the left hand side of (17), we get

 (unk+1H,1−unkH,1+unkH,2−unk−1H,2,unk+1H,1−unkH,1) ≥∥unk+1H,1−unkH,1∥2−γ∥unkH,2−unk−1H,2∥∥unk+1H,1−unkH,1∥ ≥12∥unk+1H,1−unkH,1∥2−γ22∥unkH,2−unk−1H,2∥2

For the left hand side of (18), we have

 nk+1−1∑n=nk(m(un+1H,2−unH,2)+unkH,1−unk−1H,1,un+1H,2−unH,2) ≥nk+1−1∑n=nk∥m(un+1H,2−unH,2)∥2−γ22∥unk+1H,1−unkH,1∥2−12∥unkH,2−unk−1H,2∥2 ≥m2nk+1−1∑n=nk∥un+1H,2−unH,2∥2−γ22∥unk+1H,1−unkH,1∥2

since .

Sum up the right hand side of (17) and (18), we have

 −ΔTa(unk+1H,1+unkH,2,unk+1H,1−unkH,1)−(1−ω)ΔTa(unkH,1,unk+1H,2−unkH,2) (19) −ωΔTa(unk+1H,1,unk+1H,2−unkH,2)−ΔTnk+1−1∑n=nk(unH,2,un+1H,2−unH,2) = −ΔTa(unk+1H,1,unk+1H,1−unkH,1)+ΔTa(unkH,2,unkH,1)−ΔTa(unk+1H,1,unk+1H,2) +(1−ω)ΔTa(unk+1H,1−unkH,1,unk+1H,2−unkH,2)−ΔTnk+1−1∑n=nk(unH,2,un+1H,2−unH,2) =:RHS

Note that for the terms in RHS in the above inequalities, we have

 −a(unk+1H,1,unk+1H,1−unkH,1)≤−12(∥unk+1H,1∥2a+∥unk+1H,1−unkH,1∥2a−∥unkH,1∥2a), nk+1−1∑n=nk(unH,2,un+1H,2−unH,2)≤−12nk+1−1∑n=nk(∥unH,2∥2a+∥un+1H,2−unH,2∥2a−∥un+1H,2∥2a) =−12⎛⎝∥unkH,2∥2a−∥unk+1H,2∥2a+nk+1−1∑n=nk∥un+1H,2−unH,2∥2a⎞⎠, a(unk+1H,1−unkH,1,unk+1H,2−unkH,2)≤12(∥unk+1H,1−unkH,1∥2a+∥unk+1H,2−unkH,2∥2a).

Substitute these into the left of (19) and regroup terms, we get

 RHS ≤−ΔT2∥unk+1H∥2a+ΔT2∥unkH∥2a+ΔT2nk+1−1∑n=nk∥un+1H,2−unH,2∥2a −ωΔT2∥u