Contrast-independent partially explicit time discretizations for multiscale wave problems
Eric T. Chung, Yalchin Efendiev, Wing Tat Leung, Petr N. Vabishchevich
aa r X i v : . [ m a t h . NA ] F e b Contrast-independent partially explicit time discretizations formultiscale wave problems
Eric T. Chung ∗ , Yalchin Efendiev † , Wing Tat Leung ‡ , Petr N. Vabishchevich § March 1, 2021
Abstract
In this work, we design and investigate contrast-independent partially explicit time discretizations for waveequations in heterogeneous high-contrast media. We consider multiscale problems, where the spatial het-erogeneities are at subgrid level and are not resolved. In our previous work [Chung, Efendiev, Leung,and Vabishchevich, Contrast-independent partially explicit time discretizations for multiscale flow problems,arXiv:2101.04863], we have introduced contrast-independent partially explicit time discretizations and ap-plied to parabolic equations. The main idea of contrast-independent partially explicit time discretization is tosplit the spatial space into two components: contrast dependent (fast) and contrast independent (slow) spacesdefined via multiscale space decomposition. Using this decomposition, our goal is further appropriately tointroduce time splitting such that the resulting scheme is stable and can guarantee contrast-independentdiscretization under some suitable (reasonable) conditions. In this paper, we propose contrast-independentpartially explicitly scheme for wave equations. The splitting requires a careful design. We prove that theproposed splitting is unconditionally stable under some suitable conditions formulated for the second space(slow). This condition requires some type of non-contrast dependent space and is easier to satisfy in the“slow” space. We present numerical results and show that the proposed methods provide results similar toimplicit methods with the time step that is independent of the contrast.
Multiscale problems for wave equations have been of interest for many applications. These include waveequations in materials, underground, and so on. In many applications, material properties have multiscalenature and high contrast, expressed as the ratio between media properties. These typically include somefeatures with very distinct properties in small portions of the domain, or small narrow strips (known aschannels). These heterogeneities bring challenges to numerical simulations as one needs to resolve thescales and the contrast. Recently, many spatial multiscale methods have been introduced to resolve spatialheterogeneities. However, because of high contrast, one needs to take small time step in wave equationswhen using explicit methods. To avoid this difficulty, we introduce partially explicit methods.Explicit methods are commonly used for wave equations due to the finite speed of propagation. Forexample, in seismology applications, staggered explicit methods ([39, 14]) are used, which conserve theenergy. Explicit methods have advantages as they provide a fast marching in time and can easily preserve ∗ Department of Mathematics, The Chinese University of Hong Kong (CUHK), Hong Kong SAR † Department of Mathematics, Texas A&M University, College Station, TX 77843, USA & North-Eastern Federal University,Yakutsk, Russia ‡ Department of Mathematics, University of California, Irvine, USA § Nuclear Safety Institute, Russian Academy of Sciences, Moscow, Russia & North-Eastern Federal University, Yakutsk,Russia dt = H ) and thus much coarser.We present several numerical results. We consider different heterogeneous media and different sourcefrequencies. Examples are selected where additional basis functions provide an improvement by choosing“singular” source terms, which are common for wave equations. We compare various methods and show thatthe proposed methods provide an approximation that is independent of the contrast.The paper is organized as follows. In the next section, we present Preliminaries. In Section 3, we presenta general construction of partially explicit methods. Section 4 is devoted to the construction of multiscalespaces. We present numerical results in Section 5. The conclusions are presented in Section 6. We will consider the second order wave equation in heterogeneous domain. The problem consists of finding u such that ∂ ∂t u = ∇ · ( κ ∇ u ) in Ω , (2.1)where κ ∈ L ∞ (Ω) is a high contrast heterogeneous field. The equation (2.1) is equipped with initial andboundary conditions, u (0 , · ) = u ( x ) , u t (0 , · ) = u ( x ) , and u ( t, x ) = g ( x ) on ∂ Ω .We can write the problem in the weak formulation. Find u ( t, · ) ∈ V := H (Ω) such that ( ∂ ∂t u, v ) = − a ( u, v ) ∀ v ∈ V, (2.2)where a ( u, v ) = Z Ω κ ∇ u · ∇ v. We take homogeneous boundary conditions.For semidiscretization in space, we seek an approximation u H ( t ) ∈ V H , where V H ( V H ⊂ V ) is a finitedimensional space ( H is a spatial mesh size), d dt ( u H ( t ) , v ) + a ( u H , v ) = 0 ∀ v ∈ V H , < t ≤ T, (2.3) u (0) = u H , dudt (0) = u H , (2.4)where ( u H , v ) = ( u , v ) , ( u H , v ) = ( u , v ) , ∀ ∈ V H . We set v = du H /dt in (2.3) and obtain the equalitiesfor conservation and stability for (2.3), (2.4) with respect to initial conditions: E ( t ) = E (0) , < t ≤ T, (2.5)where E ( t ) = (cid:13)(cid:13)(cid:13)(cid:13) dudt ( t ) (cid:13)(cid:13)(cid:13)(cid:13) + k u ( t ) k a , k u k a = a ( u, u ) and k u k = ( u, u ) .We take a uniform mesh with the size τ and let t n = nτ, n = 0 , . . . , N, N τ = T , u nH = u ( t n ) . Stabilityconditions for three-layer schemes with second-order accuracy with respect to τ is well known (see, forexample, [35, 36]).When using implicit discretization, u nH ≈ u H ( t n ) is sought as (cid:18) u n +1 H − u nH + u n − H τ , v (cid:19) + 14 a ( u n +1 H + 2 u nH + u n − H , v ) = 0 ∀ v ∈ V H ,n = 1 , . . . , N − , (2.6)with appropriate initial conditions. We introduce new variables s n = u nH + u n − H , r n = u nH − u n − H τ , and from (2.6), we get (cid:18) r n +1 − r n τ , v (cid:19) + 12 a ( s n +1 + s n , v ) = 0 . We take v = 2( s n +1 − s n ) = τ ( r n +1 + r n ) , which gives k r n +1 k − k r n k + k s n +1 k a − k s n k a = 0 . Thus, we have E n +1 / = E n − / , (2.7)where E n +1 / = (cid:13)(cid:13)(cid:13)(cid:13) u n +1 H − u nH τ (cid:13)(cid:13)(cid:13)(cid:13) + (cid:13)(cid:13)(cid:13)(cid:13) u n +1 H + u nH (cid:13)(cid:13)(cid:13)(cid:13) a , n = 1 , . . . , N − . This equality is a discrete version of (2.5) and guarantees unconditional stability.Stability condition for explicit method (cid:18) u n +1 H − u nH + u n − H τ , v (cid:19) + a ( u nH , v ) = 0 ∀ v ∈ V H , n = 1 , . . . , N − , (2.8)is done in a similar way. Taking into account u nH = u n +1 H + 2 u nH + u n − H − τ u n +1 H − u nH + u n − H τ , in new variables defined, (2.8) can be written as k r n +1 k − τ k r n +1 k a − k r n k + τ k r n k a + k s n +1 k a − k s n k a = 0 . Thus, we have the equality of the energies, where the energy is defined as E n +1 / = (cid:13)(cid:13)(cid:13)(cid:13) u n +1 H − u nH τ (cid:13)(cid:13)(cid:13)(cid:13) − τ (cid:13)(cid:13)(cid:13)(cid:13) u n +1 H − u nH τ (cid:13)(cid:13)(cid:13)(cid:13) a + (cid:13)(cid:13)(cid:13)(cid:13) u n +1 H + u nH (cid:13)(cid:13)(cid:13)(cid:13) a , n = 1 , . . . , N − . The quantity E n +1 / defines a norm if k v k ≥ τ k v k a ∀ v ∈ V H . (2.9)4ecause of (2.9), the stability of (2.8) takes place if τ is sufficiently small τ ≤ τ .The schemes (2.6) and (2.8) are special cases of three-layer scheme with weights: (cid:18) u n +1 H − u nH + u n − H τ , v (cid:19) + a ( σu n +1 H + (1 − σ ) u nH + σu n − H , v ) = 0 ∀ v ∈ V H ,n = 1 , . . . , N − . (2.10)When σ = 0 . we have (2.6), and if σ = 0 , we have the scheme (2.8). Unconditionally stable is the scheme(2.10), if σ ≥ . , and conditionally stable if ≤ σ < . . We have considered limit cases σ = 0 . and σ = 0 . In this section, we will discuss a temporal splitting scheme We consider V H can be decomposed into twosubspaces V H, , V H, namely, V H = V H, + V H, . We seek the solution u H = u H, + u H, and { u nH, } Nn =1 ∈ V ,H , { u nH, } Nn =1 ∈ V H, , u H = u H, + u H, ,such that ( u n +1 H − u nH + u n − H , w ) + τ a ( u n +1 H, + u n − H, + 2 u nH, , w ) = 0 ∀ w ∈ V ,H ( u n +1 H − u nH + u n − H , w ) + τ a ( ωu nH, + (1 − ω )( u n +1 H, + u n − H, ) / u nH, , w ) = 0 ∀ w ∈ V ,H . (3.1)We will consider the case ω = 1 mostly as the case ω = 0 performs similarly and is more difficult to showstability (see Appendix A). The case ω = 1 has the following form ( u n +1 H − u nH + u n − H , w ) + τ a ( u n +1 H, + u n − H, + 2 u nH, , w ) = 0 ∀ w ∈ V ,H , (3.2) ( u n +1 H − u nH + u n − H , w ) + τ a ( u nH, + u nH, , w ) = 0 ∀ w ∈ V ,H . (3.3)This is a special case of three-layer schemes, which we will investigate for stability and show that one canuse contrast independent time step.We denote the discrete energy E n + of u H is defined as E n + = k u n +1 H − u nH k + τ X i =1 , (cid:16) k u n +1 H,i k a + k u nH,i k a (cid:17) + τ a ( u n +1 H, , u nH, ) + τ a ( u n +1 H, , u nH, ) − τ k u n +1 H, − u nH, k a (3.4)or E n + = k u n +1 H − u nH k + τ (cid:16) k u n +1 H, + u nH, k a + k u nH, + u n +1 H, k a (cid:17) − τ k u n +1 H, − u nH, k a Lemma 3.1.
For u H, and u H, satisfying (3.2) and (3.3), we have E n + = E n − . Proof.
We have ( u n +1 H − u nH + u n − H , w ) + τ a ( u n +1 H, + u n − H, + 2 u nH, , w ) = 0 ∀ w ∈ V ,H , (3.5) ( u n +1 H − u nH + u n − H , w ) + τ a ( u nH, + u nH, , w ) = 0 ∀ w ∈ V ,H . (3.6)5e consider w = u n +1 H, − u n − H, in the first equation and w = u n +1 H, − u n − H, in the second equation and obtainthe following equations ( u n +1 H − u nH + u n − H , u n +1 H, − u n − H, ) + τ a ( u n +1 H, + u n − H, + 2 u nH, , u n +1 H, − u n − H, ) = 0( u n +1 H − u nH + u n − H , u n +1 H, − u n − H, ) + τ a ( u nH, + u nH, , u n +1 H, − u n − H, ) = 0 (3.7)The sum of the left hand sides can be estimated in the following way X i =1 , ( u n +1 H − u nH + u n − H , u n +1 H,i − u n − H,i ) = ( u n +1 H − u nH + u n − H , u n +1 H − u n − H ) = k u n +1 H − u nH k − k u nH − u n − H k . Next, we will estimate the right hand side of the equations. We have a ( u n +1 H, + u n − H, + 2 u nH, , u n +1 H, − u n − H, ) = 12 a ( u n +1 H, + u n − H, , u n +1 H, − u n − H, )+ a ( u nH, , u n +1 H, − u n − H, ) (3.8)and a ( u nH, + u nH, , u n +1 H, − u n − H, ) = a ( u nH, , u n +1 H, − u n − H, ) + a ( u nH, , u n +1 H, − u n − H, ) Thus, we have τ a ( u n +1 H, + u n − H, + 2 u nH, , u n +1 H, − u n − H, ) + τ a ( u nH, + u nH, , u n +1 H, − u n − H, ) = τ B + τ B + τ B where B = a ( u n +1 H, + u n − H, , u n +1 H, − u n − H, ) ,B = a ( u nH, , u n +1 H, − u n − H, ) + a ( u nH, , u n +1 H, − u n − H, ) ,B = a ( u nH, , u n +1 H, − u n − H, ) . To estimate B , we have B = a ( u n +1 H, + u n − H, , u n +1 H, − u n − H, ) = k u n +1 H, k a − k u n − H, k a =( k u n +1 H, k a + k u nH, k a ) − ( k u nH, k a + k u n − H, k a ) . We next estimate B and have B = a ( u nH, , u n +1 H, − u n − H, ) + a ( u nH, , u n +1 H, − u n − H, ) = a ( u nH, , u n +1 H, ) + a ( u nH, , u n +1 H, ) − a ( u nH, , u n − H, ) − a ( u nH, , u n − H, ) . Since a ( · , · ) is symmetric, we have B = (cid:16) a ( u n +1 H, , u nH, ) + a ( u n +1 H, , u nH, ) (cid:17) − (cid:16) a ( u nH, , u n − H, ) + a ( u nH, , u n − H, ) (cid:17) . To estimate B , we have B = a ( u nH, , u n +1 H, − u n − H, ) = a ( u nH, , u n +1 H, ) − a ( u nH, , u n − H, ) = a ( u n +1 H, , u nH, ) − a ( u nH, , u n − H, ) .
6e also have a ( u n +1 H, , u nH, ) = 12 (cid:16) k u n +1 H, k a + k u nH, k a − k u n +1 H, − u nH, k a (cid:17) and a ( u nH, , u n − H, ) = 12 (cid:16) k u nH, k a + k u n − H, k a − k u nH, − u n − H, k a (cid:17) . Thus, we have τ a ( u n +1 H, + u n − H, + 2 u nH, , u n +1 H, − u n − H, ) + τ a ( u nH, + u nH, , u n +1 H, − u n − H, ) = τ (cid:16) k u n +1 H, k a + k u nH, k a (cid:17) − τ (cid:16) k u nH, k a + k u n − H, k a (cid:17) + τ (cid:16) a ( u n +1 H, , u nH, ) + a ( u n +1 H, , u nH, ) (cid:17) − τ (cid:16) a ( u nH, , u n − H, ) + a ( u nH, , u n − H, ) (cid:17) + τ (cid:16) k u n +1 H, k a + k u nH, k a − k u n +1 H, − u nH, k a (cid:17) − τ (cid:16) k u nH, k a + k u n − H, k a − k u nH, − u n − H, k a (cid:17) (3.9)By definition of E n + , we have E n + = k u n +1 H − u nH k + X i =1 , (cid:16) τ k u n +1 H,i k a + k u nH,i k a ) (cid:17) + τ a ( u n +1 H, , u nH, ) + a ( u n +1 H, , u nH, ) − τ k u n +1 H, − u nH, k a = k u nH − u n − H k L + X i =1 , (cid:16) τ k u nH,i k a + k u n − H,i k a ) (cid:17) + τ a ( u nH, , u n − H, ) + a ( u nH, , u n − H, ) − τ k u nH, − u n − H, k a = E n − . (3.10)We will discuss two cases. In the first case, we assume V H, and V H, are orthogonal and in the secondcase, we will consider the case when they are not orthogonal. V H, and V H, are orthogonal Theorem 3.1.
The partially explicit scheme (3.2) and (3.3) is stable if k v k ≥ τ k v k a ∀ v ∈ V ,H . (3.11) Proof.
In this case, we can show that E n + = k u n +1 H − u nH k + X i =1 , (cid:16) τ k u n +1 H,i k a + k u nH,i k a ) (cid:17) + τ a ( u n +1 H, , u nH, ) + a ( u n +1 H, , u nH, ) − τ k u n +1 H, − u nH, k a ≥ k u n +1 H, − u nH, k L + X i =1 , (cid:16) τ k u n +1 H,i k a + k u nH,i k a ) (cid:17) + τ a ( u n +1 H, , u nH, ) + a ( u n +1 H, , u nH, ) (3.12)7his defines a norm. In this norm, we have the stability.We also refer to Appendix A, where we present the proof for the case ω = 0 when the spaces V H, and V H, are orthogonal. V H, and V H, are non-orthogonal We define γ, γ a < and α ∈ R as γ := sup v ∈ V H, ,v ∈ V H, ( v , v ) k v kk v k , γ a := sup v ∈ V H, ,v ∈ V H, a ( v , v ) k v k a k v k a and α = sup v ∈ V H, k v k a k v k . Lemma 3.2. If − γ ) α − ≥ τ , (3.13)we have − γ − α τ ! k u n +1 H, − u nH, k + τ (1 − γ a )2 X i =1 , (cid:16) k u n +1 H,i k a + k u nH,i k a (cid:17) ≤ E . Proof.
Since we have E = E n + for any n ≥ , we have k u n +1 H − u nH k + X i =1 , (cid:16) τ k u n +1 H,i k a + k u nH,i k a ) (cid:17) + τ a ( u n +1 H, , u nH, )+ a ( u n +1 H, , u nH, ) − τ k u n +1 H, − u nH, k a = E . (3.14)We have k u n +1 H − u nH k = X i =1 , k u n +1 H,i − u nH,i k + 2( u n +1 H, − u nH, , u n +1 H, − u nH, ) ≥ X i =1 , k u n +1 H,i − u nH,i k − γ k u n +1 H, − u nH, kk u n +1 H, − u nH, k ≥ (1 − γ ) k u n +1 H, − u nH, k . If − γ ) α − ≥ τ , we have (1 − γ ) k u n +1 H, − u nH, k − τ k u n +1 H, − u nH, k a ≥ (1 − γ − α τ k u n +1 H, − u nH, k . We also obtain that τ X i =1 , (cid:16) k u n +1 H,i k a + k u nH,i k a (cid:17) + τ a ( u n +1 H, , u nH, ) + a ( u n +1 H, , u nH, ) ≥ τ X i =1 , (cid:16) k u n +1 H,i k a + k u nH,i k a (cid:17) − γ a τ (cid:16) k u n +1 H, k a k u nH, k a + k u n +1 H, k a k u nH, k a (cid:17) γ a τ (cid:16) k u n +1 H, k a k u nH, k a + k u n +1 H, k a k u nH, k a (cid:17) ≤ γ a τ X i =1 , (cid:16) k u n +1 H,i k a + k u nH,i k a (cid:17) . Therefore, we have τ X i =1 , (cid:16) k u n +1 H,i k a + k u nH,i k a (cid:17) + τ a ( u n +1 H, , u nH, ) + a ( u n +1 H, , u nH, ) ≥ τ (1 − γ a )2 X i =1 , (cid:16) k u n +1 H,i k a + k u nH,i k a (cid:17) (3.15)and obtain E = X i =1 , (cid:16) k u n +1 H,i − u nH,i k + τ k u n +1 H,i k a + k u nH,i k a ) (cid:17) + τ a ( u n +1 H, , u nH, )+ a ( u n +1 H, , u nH, ) − τ k u n +1 H, − u nH, k a ≥ (1 − γ − α τ k u n +1 H, − u nH, k + τ (1 − γ a )2 X i =1 , (cid:16) k u n +1 H,i k a + k u nH,i k a (cid:17) . (3.16) V H, and V H, constructions In this section, we introduce a possible way to construct the spaces satisfying (3.13). Here, we follow ourprevious work [23]. We will show that the constrained energy minimization finite element space is a goodchoice of V H, since the CEM basis functions are constructed such that they are almost orthogonal to aspace ˜ V which can be easily defined. To obtain a V H, satisfying the condition (3.13), one of the possibleway is using an eigenvalue problem to construct the local basis function. Before, discussing the constructionof V H, , we will first introduce the CEM finite element space. In the following, we let V ( S ) = H ( S ) for aproper subset S ⊂ Ω . In this section, we will discuss the CEM method for solving the problem (2.2). We will construct the finiteelement space by solving a constrained energy minimization problem. We let T H be a coarse grid partitionof Ω . For each element K i ∈ T H , we consider a set of auxiliary basis functions { ψ ( i ) j } L i j =1 ∈ V ( K j ) . We thencan define a projection operator Π K i : L ( K i ) V ( i ) aux ⊂ L ( K i ) such that s i (Π i u, v ) = s i ( u, v ) ∀ v ∈ V ( i ) aux := span { ψ ( i ) j : 1 ≤ j ≤ L i } , where s i ( u, v ) = Z K i ˜ κuv. (4.1)and ˜ κ = κH − or ˜ κ = κ P i |∇ χ i | with some partition of unity χ i .We next define a projection operator by Π : L (Ω) V aux ⊂ L (Ω) s (Π u, v ) = s ( u, v ) ∀ v ∈ V aux := N e X i =1 V ( i ) , s ( u, v ) := P N e i =1 s i ( u | K i , v | K i ) . For each auxiliary basis functions ψ ( i ) j , we can define a local basisfunction φ ( i ) j ∈ V ( K + i ) such that a ( φ ( i ) j , v ) + s ( µ ( i ) j , v ) = 0 ∀ v ∈ V ( K + i ) s ( φ ( i ) j , ν ) = s ( ψ ( i ) j , ν ) ∀ ν ∈ V aux ( K + i ) where K + i is an oversampling domain of K i , which is a few coarse blocks larger than K i [10]. We then definethe space V cem as V cem := span { φ ( i ) j : 1 ≤ i ≤ N e , ≤ j ≤ L i } , where N e is the number of coarse elements. The CEM solution u cem is given by ( ∂ ∂t ( u cem ) , v ) = − a ( u cem , v ) ∀ v ∈ V cem . We remark that the V glo is a − orthogonal to a space ˜ V := { v ∈ V : Π( v ) = 0 } . We also know that V cem is closed to V glo and therefore it is almost orthogonal to ˜ V . Thus, we can choice V cem to be V H, andconstruct a space V H, in ˜ V . V H, We discuss two choices for the space V H, ⊂ ˜ V . We will present the stability properties numerically inSection 5. We will define basis functions for each coarse neighborhood ω i , which is the union of all coarse elementshaving the i -th coarse grid node. For each coarse neighborhood ω i , we consider the following eigenvalueproblem: find ( ξ ( i ) j , γ ( i ) j ) ∈ ( V ( ω i ) ∩ ˜ V ) × R , Z ω i κ ∇ ξ ( i ) j · ∇ v = γ ( i ) j H Z ω i ξ ( i ) j v, ∀ v ∈ V ( ω i ) ∩ ˜ V . (4.2)We arrange the eigenvalues by γ ( i )1 ≤ γ ( i )2 ≤ · · · . In order to obtain a reduction in error, we will select thefirst few J i dominant eigenfunctions corresponding to smallest eigenvalues of (4.2). We define V H, = span { ξ ( i ) j | ∀ ω i , ∀ ≤ j ≤ J i } . The second choice of V H, is based on the CEM type finite element space. For each coarse element K i , wewill solve an eigenvalue problem to obtain the auxiliary basis. More precisely, we find eigenpairs ( ξ ( i ) j , γ ( i ) j ) ∈ ( V ( K i ) ∩ ˜ V ) × R by solving Z K i κ ∇ ξ ( i ) j · ∇ v = γ ( i ) j Z K i ξ ( i ) j v, ∀ v ∈ V ( K i ) ∩ ˜ V . (4.3)For each K i , we choose the first few J i eigenfunctions corresponding to the smallest J i eigenvalues. Thespan of these functions form a space which is called V aux, . For each auxiliary basis function ξ ( i ) j ∈ V aux, ,10e define a basis function ζ ( i ) j ∈ V ( K + i ) such that µ ( i ) j ∈ V aux, , µ ( i ) , j ∈ V aux, and a ( ζ ( i ) j , v ) + s ( µ ( i ) , j , v ) + ( µ ( i ) , j , v ) = 0 , ∀ v ∈ V ( K + i ) , (4.4) s ( ζ ( i ) j , ν ) = 0 , ∀ ν ∈ V aux, , (4.5) ( ζ ( i ) j , ν ) = ( ξ ( i ) j , ν ) , ∀ ν ∈ V aux, . (4.6)where we use the notation V aux, to denote the space V aux defined in Section 4.1, and K + i is an oversamplingdomain a few coarse blocks larger than K i (see [10]). We define V H, = span { ζ ( i ) j | ∀ K i , ∀ ≤ j ≤ J i } . In this section, we will present representative numerical results. We consider the following mesh and timestep parameters in all examples. H = 1 / , h = 1 / , dt = 0 . , T = 0 . . We consider two medium parameter κ ( x ) , where one is a simpler compared to the other (see Figure 5.1). Bothmedium parameters are high contrast and multiscale. We choose the source term as a source distrubuted ina small region as shown in Figure 5.1). It, f , is given by f ( t, x ) = 2 − /f h exp ( − π f ( t − /f ) ) f x ( x ) , (5.1)where f is a frequency. We will consider two values for f . In numerical results, we will refer to the case(Case 1 refers to first medium and Case 2 refers to the second medium) and the frequency f . Figure 5.1: Left: κ for Case 1 (denote κ ( x ) ); Middle: κ for Case 2 (denote κ ( x ) ); Right: f x ( x ) from (5.1). f = 1 / . In our first numerical example, we consider the case with the conductivity κ ( x ) . In Figure 5.2, we plotthe reference solution computed on the fine grid (top-left), the solution computed using CEM-GMsFEM(top-right) (without additional basis functions), the solution computed using additional basis functions withimplicit method (bottom-left), and the solution computed using additional basis functions using partiallyexplicit method (bottom-right). In all cases, the reference solution is computed using standard explicit finiteelement discretization with small time step (in our case, τ = 1 e − ). First, we note that it can be seenthat additional basis functions provide improved results. This is because of the choice of source term sinceCEM-GMsFEM solution requires more detailed information to improve the solution. In Figure 5.2, we plotthe solution at the time T = 0 . . In Figure 5.3, we plot the solution and its approximations (as in Figure11.2) for T = 0 . . The errors are plotted in Figure 5.4 and Figure 5.5, where we plot both L and energyerrors. The errors are computed in a standard way using finite element discretization. In Figure 5.5, wezoom the error graph into small time interval since the error at initial time is larger. Our main observationis that our proposed approach that treats additional degrees of freedom explicitly provides a similar resultas the approach where all degrees of freedom are treated implicitly. -1.2-1-0.8-0.6-0.4-0.200.20.40.60.8 -1.2-1-0.8-0.6-0.4-0.200.20.40.60.8 -1.2-1-0.8-0.6-0.4-0.200.20.40.60.8 -1.2-1-0.8-0.6-0.4-0.200.20.40.60.8 Figure 5.2: Snapshot at T = 0 . . Top-left: reference solution. Top-Right: Implicit CEM-GMsFEM solution.Top-left: Proposed splitting method with additional basis functions. Top-right: Implicit CEM-GMsFEMwith additional basis. 12 -0.8-0.6-0.4-0.200.20.40.60.811.2 -0.8-0.6-0.4-0.200.20.40.60.811.2 -0.8-0.6-0.4-0.200.20.40.60.811.2 -0.8-0.6-0.4-0.200.20.40.60.811.2 Figure 5.3: Snapshot at T = 0 . . op-left: reference solution. Top-Right: Implicit CEM-GMsFEM solution.Top-left: Proposed splitting method with additional basis functions. Top-right: Implicit CEM-GMsFEMwith additional basis. Implicit CEM (backward Euler)Implicit CEM with Additional BasisPartially Explicit Splitting CEM
Implicit CEM (backward Euler)Implicit CEM with Additional BasisPartially Explicit Splitting CEM
Figure 5.4: Second type of V ,H (CEM Dof: , V ,H Dof: ). Left: L error. Right: Energy error.13 .2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.600.050.10.150.20.250.3 Implicit CEM (backward Euler)Implicit CEM with Additional BasisPartially Explicit Splitting CEM
Implicit CEM (backward Euler)Implicit CEM with Additional BasisPartially Explicit Splitting CEM
Figure 5.5: The error in the time interval from . to . . Second type of V ,H is used (CEM Dof: , V ,H Dof: ). Left: L error. Right: Energy error. f = 1 / . In our second numerical example, we consider the case with the conductivity κ ( x ) . As in the previous case,we first depict the reference solution (top-left), the solution computed using CEM-GMsFEM (top-right)(without additional basis functions), the solution computed using additional basis functions with implicitmethod (bottom-left), and the solution computed using additional basis functions using partially explicitmethod (bottom-right) in Figure 5.6. We see that additional basis functions provide an improved result.Moreover, there is a little difference between the implicit solution that uses additional basis functions andthe partially explicit solution that treats additional degrees of freedom explicitly. In Figure 5.6, we plot thesolution at the time T = 0 . and we can make similar observations. In Figure 5.7, we plot the solution andits approximations (as in Figure 5.6) for T = 0 . . The errors are plotted in Figure 5.8 and Figure 5.9, wherewe plot both L and energy errors. In Figure 5.9, we zoom the error graph into smaller time interval. Ourmain observation is that our proposed approach that treats additional degrees of freedom explicitly providesa similar result as the approach where all degrees of freedom are treated implicitly.14 -1.5-1-0.500.511.5 -0.500.511.5 -0.500.511.5 -0.500.511.5 Figure 5.6: Snapshot at T = 0 . . Top-left: reference solution. Top-Right: Implicit CEM-GMsFEM solution.Top-left: Proposed splitting method with additional basis functions. Top-right: Implicit CEM-GMsFEMwith additional basis. -1.5-1-0.500.511.5 -1.5-1-0.500.511.5 -1.5-1-0.500.511.5 -1.5-1-0.500.511.5 Figure 5.7: Snapshot at T = 0 . . op-left: reference solution. Top-Right: Implicit CEM-GMsFEM solution.Top-left: Proposed splitting method with additional basis functions. Top-right: Implicit CEM-GMsFEMwith additional basis. 15 Implicit CEM (backward Euler)Implicit CEM with Additional BasisPartially Explicit Splitting CEM
Implicit CEM (backward Euler)Implicit CEM with Additional BasisPartially Explicit Splitting CEM
Figure 5.8: Second type of V ,H (CEM Dof: , V ,H Dof: ). Left: L error. Right: Energy error. Implicit CEM (backward Euler)Implicit CEM with Additional BasisPartially Explicit Splitting CEM
Implicit CEM (backward Euler)Implicit CEM with Additional BasisPartially Explicit Splitting CEM
Figure 5.9: The error in the time interval from . to . . Second type of V ,H is used (CEM Dof: , V ,H Dof: ). Left: L error. Right: Energy error. f = 1 In our final numerical example, we consider the case with the conductivity κ ( x ) and higher frequency. Inthis case, the error is larger as expected and additional basis functions in CEM-GMsFEM provide improvedsolution. As in the previous cases, we first depict the reference solution (top-left), the solution computed usingCEM-GMsFEM (top-right) (without additional basis functions), the solution computed using additional basisfunctions with implicit method (bottom-left), and the solution computed using additional basis functionsusing partially explicit method (bottom-right) in Figure 5.10. We observe that there is a little differencebetween the implicit solution that uses additional basis functions and the partially explicit solution thattreats additional degrees of freedom explicitly. In Figure 5.10, we plot the solution at the time T = 0 . andwe can make similar observations. In Figure 5.11, we plot the solution and its approximations (as in Figure5.10) for T = 0 . . The errors are plotted in Figure 5.12 and Figure 5.13, where we plot both L and energyerrors. In Figure 5.13, we zoom the error graph into smaller time interval. Our main observation is thatour proposed approach that treats additional degrees of freedom explicitly provides a similar result as theapproach where all degrees of freedom are treated implicitly.16 -0.4-0.3-0.2-0.100.10.20.30.4 -0.4-0.3-0.2-0.100.10.20.30.4 -0.4-0.3-0.2-0.100.10.20.30.4 -0.4-0.3-0.2-0.100.10.20.30.4 Figure 5.10: Snapshot at T = 0 . . Top-left: reference solution. Top-Right: Implicit CEM-GMsFEM solution.Top-left: Proposed splitting method with additional basis functions. Top-right: Implicit CEM-GMsFEMwith additional basis. -0.4-0.3-0.2-0.100.10.20.3 -0.4-0.3-0.2-0.100.10.20.3 -0.4-0.3-0.2-0.100.10.20.3 -0.4-0.3-0.2-0.100.10.20.3 Figure 5.11: Snapshot at T = 0 . . op-left: reference solution. Top-Right: Implicit CEM-GMsFEM solution.Top-left: Proposed splitting method with additional basis functions. Top-right: Implicit CEM-GMsFEMwith additional basis. 17 Implicit CEM (backward Euler)Implicit CEM with Additional BasisPartially Explicit Splitting CEM
Implicit CEM (backward Euler)Implicit CEM with Additional BasisPartially Explicit Splitting CEM
Figure 5.12: Second type of V ,H (CEM Dof: , V ,H Dof: ). Left: L error. Right: Energy error. Implicit CEM (backward Euler)Implicit CEM with Additional BasisPartially Explicit Splitting CEM
Implicit CEM (backward Euler)Implicit CEM with Additional BasisPartially Explicit Splitting CEM
Figure 5.13: The error in the time interval from . to . . Second type of V ,H is used (CEM Dof: , V ,H Dof: ). Left: L error. Right: Energy error. The proposed methods are still coupled via mass matrix. One way to remove this coupling is via masslumping methods [15]. There are many ways to do mass lumping, which can be applied in this paper. Wehave tried some mass lumping schemes; however, their performance were not optimal. We believe for masslumping other discretizations (mixed or Discontinuous Galerkin, e.g., [5]) may be more effective. Here, weconsider one example with mass lumping.In this subsection, we will introduce a way for mass laamping. We will consider the auxiliary basisfunctions of V H, are defined as following: for each coarse element K i ∈ T H , we consider ψ ( i ) j = I K ( i ) j ( x ) where I K ( i ) j is the characteristic function of K ( i ) j ⊂ K i . K ( i ) j is defined as K ( i )1 is the region with low wavespeed and K ( i )2 is the region with high wave speed. For example, we consider K ( i )1 = { x ∈ K i | κ ( x ) ≤ } K ( i )2 = { x ∈ K i | κ ( x ) > } . For each coarse element K i , we can solve the eigenvalue problem to obtain the auxiliary basis for V ,H . Wecan find ( ξ ( i ) j , γ ( i ) j ) ∈ ( V ( K i ) ∩ ˜ V ) × R , Z ω i κ ∇ ξ ( i ) j · ∇ v = γ ( i ) j Z ω i ξ ( i ) j v ∀ v ∈ V ( K i ) ∩ ˜ V ˜ V = { v ∈ V | ( v, ψ ( i ) j ) = 0 ∀ i, j } . The auxiliary space V aux, and V aux, are defined as V aux, = span i,j { ψ ( i ) j } ,V aux, = span i,j { ξ ( i ) j } . The multiscale basis functions φ ( i ) j, of V H, are obtained by finding ( φ ( i ) j, , µ ( i ) j, ) ∈ V ( K + i ) × ( V aux, + V aux, ) a ( φ ( i ) j, , v ) + ( µ ( i ) j, , v ) = 0 ∀ v ∈ V ( K + i )( φ ( i ) j, , ψ ( l ) k ) = δ il δ jk ( φ ( i ) j, , ξ ( l ) k ) = 0 . Similarly, the multiscale basis functions φ ( i ) j, of V H, are obtained by finding ( φ ( i ) j, , µ ( i ) j, ) ∈ V ( K + i ) × ( V aux, + V aux, ) a ( φ ( i ) j, , v ) + ( µ ( i ) j, , v ) = 0 ∀ v ∈ V ( K + i )( φ ( i ) j, , ψ ( l ) k ) = 0( φ ( i ) j, , ξ ( l ) k ) = δ il δ jk . The multiscale finite element spaces V H, and V H, are defined as V H, = span i,j { φ ( i ) j, } ,V H, = span i,j { φ ( i ) j, } . For u H ∈ V H, + V H, , we have u H = P i,j u ( i ) j, φ ( i ) j, + P i,j u ( i ) j, φ ( i ) j, and a ( u, v ) = − ( X i,j u ( i ) j, µ ( i ) j, , v ) − ( X i,j u ( i ) j, µ ( i ) j, , v ) ∀ v ∈ V H, + V H, . Thus, we can consider the weak formulation of u H,tt = ∇ · ( κ ∇ u H ) + f as ( X i,j u ( i ) j, φ ( i ) j, + X i,j u ( i ) j, φ ( i ) j, , w ) = ( X i,j u ( i ) j, µ ( i ) j, + X i,j u ( i ) j, µ ( i ) j, + f, w ) ∀ w ∈ W H for some testing space W H . We can consider the W H to be ( V aux, + V aux, ) and we have ( X i,j ( u ( i ) j, ) tt φ ( i ) j, + X i,j ( u ( i ) j, ) tt φ ( i ) j, , ψ ( l ) k ) = ( X i,j u ( i ) j, µ ( i ) j, + X i,j u ( i ) j, µ ( i ) j, + f, ψ ( l ) k ) ∀ l, k ( X i,j ( u ( i ) j, ) tt φ ( i ) j, + X i,j ( u ( i ) j, ) tt φ ( i ) j, , ξ ( l ) k ) = ( X i,j u ( i ) j, µ ( i ) j, + X i,j u ( i ) j, µ ( i ) j, + f, ξ ( l ) k ) ∀ l, k. Since ( φ ( i ) j, , ψ ( l ) k ) = 0( φ ( i ) j, , ξ ( l ) k ) = δ il δ jk , and ( φ ( i ) j, , ψ ( l ) k ) = δ il δ jk ( φ ( i ) j, , ξ ( l ) k ) = 0 ,
19e have ( u ( l ) k, ) tt = ( X i,j u ( i ) j, µ ( i ) j, + X i,j u ( i ) j, µ ( i ) j, + f, ψ ( l ) k ) ∀ l, k ( u ( l ) k, ) tt = ( X i,j u ( i ) j, µ ( i ) j, + X i,j u ( i ) j, µ ( i ) j, + f, ξ ( l ) k ) ∀ l, k. Since ( X i,j u ( i ) j, µ ( i ) j, + X i,j u ( i ) j, µ ( i ) j, + f, v ) = − a ( u H , v ) , and ( X i,j u ( i ) j, µ ( i ) j, + X i,j u ( i ) j, µ ( i ) j, , ψ ( l ) k ) = ( X i,j u ( i ) j, µ ( i ) j, + X i,j u ( i ) j, µ ( i ) j, , φ ( l ) k, ) = − a ( u H , φ ( l ) k, ) , ( X i,j u ( i ) j, µ ( i ) j, + X i,j u ( i ) j, µ ( i ) j, , ξ ( l ) k ) = ( X i,j u ( i ) j, µ ( i ) j, + X i,j u ( i ) j, µ ( i ) j, , φ ( l ) k, ) = − a ( u H , φ ( l ) k, ) , We have ( u ( l ) k, ) tt + X i,j u ( i ) j, a ( φ ( i ) j, , φ ( l ) k, ) + X i,j u ( i ) j, a ( φ ( i ) j, , φ ( l ) k, ) = ( f, ψ ( l ) k ) ∀ l, k ( u ( l ) k, ) tt + X i,j u ( i ) j, a ( φ ( i ) j, , φ ( l ) k, ) + X i,j u ( i ) j, a ( φ ( i ) j, , φ ( l ) k, ) = ( f, ξ ( l ) k ) ∀ l, k. Using discretization, we can obtain a diagonal mass matrix. We remark that this can be considered as a masslumping trick for the system (3.1). Above, we derived this for global basis functions; however, the calculationscan be localized. Instead of using the L inner product for the mass matrix, we can use ( π ( u ) , π ( v )) to obtainan approximated mass matrix where π is the L projection operator from V H to V aux, + V aux, . We haveused this localized projection idea in our numerical simulations.Following the idea of the CEM-GMsFEM, insteading of using solving a global problem to obtain the globalbasis functions, we can solve the problem in an oversampling domain to obtain the CEM basis functions.Since the CEM basis functions are exponentially converging to the global basis functions, we can use theinner product ( π ( u ) , π ( v )) for mass lumping of this CEM space.We will show a numerical example for example case 2 with f = 1 / . without mass lumpingwith mass lumping without mass lumpingwith mass lumping Figure 5.14: Mass lumping result for third type of space and τ = 0 . . (CEM Dof: , V ,H Dof: ).Left: L error , Right: Energy error 20 .2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.600.050.10.150.20.250.30.35 without mass lumpingwith mass lumping without mass lumpingwith mass lumping Figure 5.15: Mass lumping result for third type of space and τ = 0 . . Time from . to . . (CEM Dof: , V ,H Dof: ). Left: L error , Right: Energy error In this paper, we design contrast-independent partially explicit time discretization methods for wave equa-tions. The proposed methods differ from our previous works, where we first introduced contrast-independentpartial explicit methods for parabolic equations [23]. The proposed approach uses temporal splitting basedon spatial multiscale splitting. We first introduce two spatial spaces, first account for spatial features relatedto fast time scales and the second for spatial features related to “slow” time scales. Using these spaces, wepropose time splitting, where the first equation solves for fast components implicitly and the second equationsolves for slow components explicitly. Our proposed method is still implicit via mass matrix; however, itis explicit in terms of stiffness matrix for the slow component (which is contrast independent). Via masslumping, one can remove the coupling, which is briefly discussed in the paper. We show a stability of theproposed splitting under suitable conditions for the second space, where the fast components are absent.We present numerical results, which show that the proposed methods provide very similar results as fullyimplicit methods using explicit methods with the time stepping that is independent of the contrast.
A Proof for ω = 0 For the case ω = 0 , we only consider V H, ⊥ V H, . The scheme for ω = 0 has the form ( u n +1 H, − u nH, + u n − H, , w ) + τ a ( u n +1 H, + u n − H, + 2 u nH, , w ) = 0 ∀ w ∈ V H, ( u n +1 H, − u nH, + u n − H, , w ) + τ a ( u n +1 H, + u n − H, + 2 u nH, , w ) = 0 ∀ w ∈ V H, . We define an inner product ( · , · ) m τ such that ( u, v ) m τ = ( u, v ) + τ a ( u, v ) . We then define two operators b τ : V H, V H, , c τ : V H, + V H, V H, such that ( b τ ( v ) , v ) m τ = ( v , v ) ∀ v ∈ V H, , ( c τ ( u ) , v ) m τ := τ a ( u, v ) ∀ v ∈ V H, ,a ( d τ ( u ) , v ) := a ( u, v ) ∀ v ∈ V H, .
21e remark that ( c τ ( v ) , v ) m τ = ( v , v ) m τ ∀ v ∈ V H, , v ∈ V H, , and k b τ ( v ) k m τ = ( v , b τ ( v )) ≤ k v kk b τ ( v ) k≤ k v kk b τ ( v ) k . Equality holds if and only if v is a constant function.We define an inner product ( · , · ) s τ such that ( v , v ) s τ = ( v , v ) − ( b τ ( v ) , b τ ( v )) m τ for v , v ∈ V H, . We then define three (semi) norms k · k m τ , k · k s τ , k · k n τ such that k v k m τ = ( v, v ) m τ , k v k s τ = ( v , v ) − k b τ ( v ) k m τ , k v k n τ = τ k v k a − k c τ ( v ) k m τ − k d τ ( v ) k s τ . It is clear that k · k m τ , k · k s τ are (semi) norms. To show that k · k n τ is a norm in V H, , we only need tocheck k v k n τ > ∀ v ∈ V H, . Since a ( d τ ( v ) , v − d τ ( v )) = 0 , we have τ k v k a = τ k v − d τ ( v ) k a + τ k d τ ( v ) k a . We then have τ k d τ ( v ) k a − k d τ ( v ) k s τ = τ k d τ ( v ) k a − ( d τ ( v ) , d τ ( v )) + ( b τ ( d τ ( v )) , b τ ( d τ ( v ))) m τ = τ k d τ ( v ) k a − ( d τ ( v ) , d τ ( v )) + ( b τ ( d τ ( v )) , d τ ( v )) . Since ( b τ ( d τ ( v )) , d τ ( v )) − ( d τ ( v ) , d τ ( v )) = − τ a ( b τ ( d τ ( v )) , d τ ( v )) + ( b τ ( d τ ( v )) , d τ ( v )) m τ − ( d τ ( v ) , d τ ( v )) = − τ a ( b τ ( d τ ( v )) , d τ ( v )) , we have τ k d τ ( v ) k a − k d τ ( v ) k s τ = τ a ( d τ ( v ) − b τ ( d τ ( v )) , d τ ( v )) . Since ( c τ ( v ) , c τ ( v )) m τ = τ a ( c τ ( v ) , v ) = τ a ( c τ ( v ) , d τ ( v )) =( c τ ( v ) , d τ ( v )) m τ − ( c τ ( v ) , d τ ( v )) =( c τ ( v ) , d τ ( v ) − b τ d τ ( v )) m τ = τ a ( v , d τ ( v ) − b τ d τ ( v )) = τ a ( d τ ( v ) − b τ ( d τ ( v )) , d τ ( v )) ,
22e have τ k d τ ( v ) k a − k d τ ( v ) k s τ − ( c τ ( v ) , c τ ( v )) m τ = 0 . Therefore, we have k v k n τ = τ k v − d τ ( v ) k a > . Lemma A.1.
For u ∈ V H, and u ∈ V H, , we have k u k + τ k u k a − k b τ ( u ) + c τ ( u ) k m τ = k u + d τ ( u ) k s τ + k u k n τ . Proof.
We have k u k + τ k u k a − k b τ ( u ) − c τ ( u ) k m τ = k u k − k b τ ( u ) k m τ + 2( b τ ( u ) , c τ ( u )) m τ + τ k u k a − k c τ ( u ) k m τ . By definition of ( · , · ) s τ , we have k u k s τ = k u k − k b τ ( u ) k m τ . ( u , d τ ( u )) s τ = ( u , d τ ( u )) − ( b τ u , b τ d τ ( u )) m τ =( u , d τ ( u )) − ( b τ u , d τ ( u )) =( u , d τ ( u )) − ( b τ u , d τ ( u )) m τ + τ a ( b τ u , d τ ( u )) . Since ( u , d τ ( u )) = ( b τ u , d τ ( u )) m τ , we have ( u , d τ ( u )) s τ = τ a ( b τ u , d τ ( u )) = τ a ( b τ u , u ) =( b τ ( u ) , c τ ( u )) m τ . Thus, we have k u k + τ k u k a − k b τ ( u ) − c τ ( u ) k m τ = k u k s τ + 2( u , d τ ( u )) s τ + τ k u k a − k c τ ( u ) k m τ = k u + d τ ( u ) k s τ + τ k u k a − k c τ ( u ) k m τ − k d τ ( u ) k s τ = k u + d τ ( u ) k s τ + k u k n τ . We then define the energy E n + by E n + := k b ( u n +1 H, ) − c ( u n +1 H, ) − b ( u nH, ) + c ( u nH, ) k m τ + k u n +1 H, − u nH, k − τ k u n +1 H, − u nH, k a + k u n +1 H, + d τ ( u n +1 H, ) k s τ + k u n +1 H, k n τ + k u nH, + d τ ( u nH, ) k s τ + k u nH, k n τ . heorem A.1. If V H, ⊥ V H, , ω = 0 , we have E n + = E n − and the scheme is stable if sup v ∈ V H, k v k a k v k ≤ τ . Proof.
We consider the test functions b τ ( u n +1 H, − u n − H, ) − c τ ( u n +1 H, − u n − H, ) and u n +1 H, − u n − H, . We have ( u n +1 H, − u nH, + u n − H, , b τ ( u n +1 H, − u n − H, ) − c τ ( u n +1 H, − u n − H, ))+ τ a ( u n +1 H, + u n − H, + 2 u nH, , b τ ( u n +1 H, − u n − H, ) − c τ ( u n +1 H, − u n − H, )) = 0( u n +1 H, − u nH, + u n − H, , u n +1 H, − u n − H, ) + τ a ( u n +1 H, + u n − H, + 2 u nH, , u n +1 H, − u n − H, ) = 0 . We consider B = ( u n +1 H, + u n − H, , b τ ( u n +1 H, − u n − H, )) + τ a ( u n +1 H, + u n − H, , b τ ( u n +1 H, − u n − H, )) ,B = τ a ( u nH, , b τ ( u n +1 H, − u n − H, )) = 2( c τ ( u nH, ) , b τ ( u n +1 H, − u n − H, )) m τ ,B = − τ a ( u nH, , c τ ( u n +1 H, − u n − H, )) = − c τ ( u nH, ) , c τ ( u n +1 H, − u n − H, )) m τ ,B = − u nH, , b τ ( u n +1 H, − u n − H, )) = − b τ ( u nH, ) , b τ ( u n +1 H, − u n − H, )) m τ ,B = 2( u nH, , c τ ( u n +1 H, − u n − H, )) = 2( b τ ( u nH, ) , c τ ( u n +1 H, − u n − H, )) m τ ,B = − ( u n +1 H, + u n − H, , c τ ( u n +1 H, − u n − H, )) ,B = − τ a ( u n +1 H, + u n − H, , c τ ( u n +1 H, − u n − H, )) ,C = ( u n +1 H, − u nH, + u n − H, , u n +1 H, − u n − H, ) + τ a ( u nH, , u n +1 H, − u n − H, ) ,C = τ a ( u n +1 H, + u n − H, , u n +1 H, − u n − H, ) , such that X i =1 B i = X i =1 C i = 0 . By the definition of B , B , B , B , we have B + B + B + B = − b τ ( u n +1 H, ) − c τ ( u n +1 H, ) , b τ ( u nH, ) − c τ ( u nH, )) m τ . Since ( u, v ) m τ = ( u, v ) + τ a ( u, v ) , we have B =( u n +1 H, + u n − H, , b τ ( u n +1 H, − u n − H, )) + τ a ( u n +1 H, + u n − H, , b τ ( u n +1 H, − u n − H, )) =( u n +1 H, + u n − H, , b τ ( u n +1 H, − u n − H, )) m τ = ( u n +1 H, + u n − H, , u n +1 H, − u n − H, ) = k u n +1 H, k + k u nH, k − k u nH, k + k u n − H, k . =( u n +1 H, − u nH, + u n − H, , u n +1 H, − u n − H, ) + τ a ( u nH, , u n +1 H, − u n − H, ) = k u n +1 H, − u nH, k + τ k u n +1 H, k a + τ k u nH, k a − τ k u n +1 H, − u nH, k a − (cid:16) k u nH, − u n − H, k + τ k u nH, k a + τ k u n − H, k a − τ k u nH, − u n − H, k a (cid:17) .C = τ a ( u n +1 H, + u n − H, , u n +1 H, − u n − H, ) = ( u n +1 H, + u n − H, , c τ ( u n +1 H, − u n − H, )) m τ =( u n +1 H, + u n − H, , c τ ( u n +1 H, − u n − H, )) + τ a ( u n +1 H, + u n − H, , c τ ( u n +1 H, − u n − H, )) = − B − B . Therefore, we have k u n +1 H, k + k u nH, k + k u n +1 H, − u nH, k + τ k u n +1 H, k a + τ k u nH, k a − τ k u n +1 H, − u nH, k a − b τ ( u n +1 H, ) − c τ ( u n +1 H, ) , b τ ( u nH, ) − c τ ( u nH, )) m τ = k u nH, k + k u n − H, k + k u nH, − u n − H, k + τ k u nH, k a + τ k u n − H, k a − τ k u nH, − u n − H, k a − b τ ( u nH, ) − c τ ( u nH, ) , b τ ( u n − H, ) − c τ ( u n − H, )) m τ . We can observe that − b τ ( u n +1 H, ) − c τ ( u n +1 H, ) , b τ ( u nH, ) − c τ ( u nH, )) m τ = k b τ ( u n +1 H, ) − c τ ( u n +1 H, ) − b τ ( u nH, ) + c τ ( u nH, ) k m τ − k b τ ( u n +1 H, ) − c τ ( u n +1 H, ) k m τ − k b τ ( u nH, ) − c τ ( u nH, ) k m τ Thus, we have E n + = E n − . Next, we will do some formal calculations to show that our proposed energy is close to the continuousenergy when τ is small. We remark that for τ → , we have b τ ( u ) = u + O ( τ ) in L c τ ( u ) = O ( τ ) in L , k u k s τ = ( u , u ) − k b τ ( u ) k m τ ≈ τ k u k a + O ( τ ) , k u k m τ ≈ k u k + O ( τ ) , k u k n τ = τ k u k a − k c τ ( u ) k m τ − k d τ ( u ) k s τ = τ (cid:16) k u k a − k d τ ( u ) k a (cid:17) + O ( τ ) , n + := k b ( u n +1 H, ) − c ( u n +1 H, ) − b ( u nH, ) + c ( u nH, ) k m τ + k u n +1 H, − u nH, k − τ k u n +1 H, − u nH, k a + k u n +1 H, + d τ ( u n +1 H, ) k s τ + k u n +1 H, k n τ + k u nH, + d τ ( u nH, ) k s τ + k u nH, k n τ . If k u n +1 α − u nα k a = k u n +1 α − u nα k = O ( τ ) , we have E n + = k u n +1 H, − u nH, k + k u n +1 H, − u nH, k + τ k u n +1 H, + d τ ( u n +1 H, ) k a + τ (cid:16) k u n +1 H, k a − k d τ ( u n +1 H, ) k a (cid:17) + τ k u nH, + d τ ( u nH, ) k a + τ (cid:16) k u nH, k a − k d τ ( u nH, ) k a (cid:17) + O ( τ ) . Since a ( u n +1 H, + u n +1 H, , u n +1 H, + u n +1 H, ) = a ( u n +1 H, , u n +1 H, ) + 2 a ( d τ ( u n +1 H, ) , u n +1 H, ) + k u n +1 H, k a = k u n +1 H, + d τ ( u n +1 H, ) k a + (cid:16) k u n +1 H, k a − k d τ ( u n +1 H, ) k a (cid:17) , we have E n + = k u n +1 H, − u nH, k + k u n +1 H, − u nH, k + τ k u n +1 H, + u n +1 H, k a + τ (cid:16) k u nH, + u nH, k a (cid:17) + O ( τ ) and τ − E n + ≈ k u n +1 H − u nH τ k + τ (cid:16) k u n +1 H k a + k u nH k a (cid:17) . The Appendix shows that our proposed approach can also be used in splitting method with ω = 0 . References [1] A. Abdulle. Explicit methods for stiff stochastic differential equations. In
Numerical Analysis of Mul-tiscale Computations , pages 1–22. Springer, 2012.[2] G. Allaire and R. Brizzi. A multiscale finite element method for numerical homogenization.
SIAM J.Multiscale Modeling and Simulation , 4(3):790–812, 2005.[3] G. Ariel, B. Engquist, and R. Tsai. A multiscale method for highly oscillatory ordinary differentialequations with resonance.
Mathematics of Computation , 78(266):929–956, 2009.[4] U. M. Ascher, S. J. Ruuth, and R. J. Spiteri. Implicit-explicit Runge-Kutta methods for time-dependentpartial differential equations.
Applied Numerical Mathematics , 25(2-3):151–167, 1997.[5] S. W. Cheung, E. T. Chung, Y. Efendiev, and W. T. Leung. Explicit and energy-conserving con-straint energy minimizing generalized multiscale discontinuous galerkin method for wave propagationin heterogeneous media. arXiv preprint arXiv:2009.00991 , 2020.[6] E. T. Chung, Y. Efendiev, and T. Hou. Adaptive multiscale model reduction with generalized multiscalefinite element methods.
Journal of Computational Physics , 320:69–95, 2016.[7] E. T. Chung, Y. Efendiev, and C. Lee. Mixed generalized multiscale finite element methods andapplications.
SIAM Multiscale Model. Simul. , 13:338–366, 2014.[8] E. T. Chung, Y. Efendiev, and W. T. Leung. Generalized multiscale finite element methods for wavepropagation in heterogeneous media.
Multiscale Modeling & Simulation , 12(4):1691–1721, 2014.[9] E. T. Chung, Y. Efendiev, and W. T. Leung. Residual-driven online generalized multiscale finite elementmethods.
Journal of Computational Physics , 302:176–190, 2015.2610] E. T. Chung, Y. Efendiev, and W. T. Leung. Constraint energy minimizing generalized multiscale finiteelement method.
Computer Methods in Applied Mechanics and Engineering , 339:298–319, 2018.[11] E. T. Chung, Y. Efendiev, and W. T. Leung. Constraint energy minimizing generalized multiscale finiteelement method in the mixed formulation.
Computational Geosciences , 22(3):677–693, 2018.[12] E. T. Chung, Y. Efendiev, and W. T. Leung. Fast online generalized multiscale finite element methodusing constraint energy minimization.
Journal of Computational Physics , 355:450–463, 2018.[13] E. T. Chung, Y. Efendiev, W. T. Leung, M. Vasilyeva, and Y. Wang. Non-local multi-continua upscalingfor flows in heterogeneous fractured media.
Journal of Computational Physics , 372:22–34, 2018.[14] E. T. Chung, C. Y. Lam, and J. Qian. A staggered discontinuous Galerkin method for the simulationof seismic waves with surface topography.
Geophysics , 80(4):T119–T135, 2015.[15] G. Cohen, P. Joly, J. E. Roberts, and N. Tordjman. Higher order triangular finite elements with masslumping for the wave equation.
SIAM Journal on Numerical Analysis , 38(6):2047–2078, 2001.[16] L. Durlofsky. Numerical calculation of equivalent grid block permeability tensors for heterogeneousporous media.
Water Resour. Res. , 27:699–708, 1991.[17] W. E and B. Engquist. Heterogeneous multiscale methods.
Comm. Math. Sci. , 1(1):87–132, 2003.[18] Y. Efendiev, J. Galvis, and T. Hou. Generalized multiscale finite element methods (GMsFEM).
Journalof Computational Physics , 251:116–135, 2013.[19] Y. Efendiev and T. Hou.
Multiscale Finite Element Methods: Theory and Applications , volume 4 of
Surveys and Tutorials in the Applied Mathematical Sciences . Springer, New York, 2009.[20] Y. Efendiev, S. Pun, and P. N. Vabishchevich. Temporal splitting algorithms for non-stationary multi-scale problems. arXiv preprint , 2020.[21] Y. Efendiev and P. N. Vabishchevich. Splitting methods for solution decomposition in nonstationaryproblems. arXiv preprint arXiv:2008.08111 , 2020.[22] B. Engquist and Y.-H. Tsai. Heterogeneous multiscale methods for stiff ordinary differential equations.
Mathematics of computation , 74(252):1707–1742, 2005.[23] W. T. L. Eric Chung, Yalchin Efendiev and P. N. Vabishchevich. Contrast-independent partially explicittime discretizations for multiscale flow problems. arXiv:2101.04863.[24] P. Henning, A. Målqvist, and D. Peterseim. A localized orthogonal decomposition method for semi-linear elliptic problems.
ESAIM: Mathematical Modelling and Numerical Analysis , 48(5):1331–1349,2014.[25] T. Hou and X. Wu. A multiscale finite element method for elliptic problems in composite materials andporous media.
J. Comput. Phys. , 134:169–189, 1997.[26] T. Y. Hou, D. Huang, K. C. Lam, and P. Zhang. An adaptive fast solver for a general class of positivedefinite matrices via energy decomposition.
Multiscale Modeling & Simulation , 16(2):615–678, 2018.[27] T. Y. Hou, Q. Li, and P. Zhang. Exploring the locally low dimensional structure in solving randomelliptic pdes.
Multiscale Modeling & Simulation , 15(2):661–695, 2017.[28] T. Y. Hou, D. Ma, and Z. Zhang. A model reduction method for multiscale elliptic pdes with randomcoefficients using an optimization approach.
Multiscale Modeling & Simulation , 17(2):826–853, 2019.[29] P. Jenny, S. Lee, and H. Tchelepi. Multi-scale finite volume method for elliptic problems in subsurfaceflow simulation.
J. Comput. Phys. , 187:47–67, 2003.2730] T. Li, A. Abdulle, et al. Effectiveness of implicit methods for stiff stochastic differential equations. In
Commun. Comput. Phys . Citeseer, 2008.[31] G. I. Marchuk. Splitting and alternating direction methods.
Handbook of numerical analysis , 1:197–462,1990.[32] H. Owhadi and L. Zhang. Metric-based upscaling.
Comm. Pure. Appl. Math. , 60:675–723, 2007.[33] A. Roberts and I. Kevrekidis. General tooth boundary conditions for equation free modeling.
SIAM J.Sci. Comput. , 29(4):1495–1510, 2007.[34] G. Samaey, I. Kevrekidis, and D. Roose. Patch dynamics with buffers for homogenization problems.
J.Comput. Phys. , 213(1):264–287, 2006.[35] A. A. Samarskii.
The Theory of Difference Schemes . Marcel Dekker, New York, 2001.[36] A. A. Samarskii, P. P. Matus, and P. N. Vabishchevich.
Difference Schemes with Operator Factors .Kluwer Academic Pub, 2002.[37] P. N. Vabishchevich. Splitting methods for dynamic second order equations. submitted.[38] P. N. Vabishchevich.
Additive Operator-Difference Schemes: Splitting Schemes . Walter de GruyterGmbH, Berlin, Boston, 2013.[39] J. Virieux. Sh-wave propagation in heterogeneous media: Velocity-stress finite-difference method.
Geo-physics , 49(11):1933–1942, 1984.[40] X. Wu, Y. Efendiev, and T. Hou. Analysis of upscaling absolute permeability.