{'fraction_non_alphanumeric': 0.043989138484324855, 'fraction_numerical': 0.011799555665267835, 'mean_word_length': 4.066533266633317, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Computational multiscale methods for first-order wave equation using mixed CEM-GMsFEM

In this paper, we consider a pressure-velocity formulation of the heterogeneous wave equation and employ the constraint energy minimizing generalized multiscale finite element method (CEM-GMsFEM) to solve this problem. The proposed method provides a flexible framework to construct crucial multiscale basis functions for approximating the pressure and velocity. These basis functions are constructed by solving a class of local auxiliary optimization problems over the eigenspaces that contain local information on the heterogeneity. Techniques of oversampling are adapted to enhance the computational performance. The first-order convergence of the proposed method is proved and illustrated by several numerical tests.


Introduction
Wave propagation and its numerical simulations have been widely studied for years due to its fundamental importance in engineering applications. For example, these problems arise in the study of seismic wave propagation from geoscience [26]. In such applications, the background materials in the domain are often highly heterogeneous, and their elastic properties may vary with the depth rapidly. In those cases with non-smooth heterogeneous coefficients, direct simulation using standard numerical methods (e.g. finite element method [20]) may lead to prohibitively expensive computational cost in order to resolve the heterogeneous structure of the media, though traditional methods can capture the fine-scale features with moderately high computational resources [12]. Therefore, it is necessary to apply model reduction techniques to alleviate the computational burden in the accurate simulations of wave propagation.
Many model reduction techniques have been well developed in the existing literature. For example, in numerical upscaling methods [18,25,27,28], one typically derives some upscaling media and solves the resulting upscaled problem globally on a coarse grid. The dimensions of the corresponding linear systems are much smaller, giving a guaranteed saving of computational cost. In addition, various multiscale methods [5,15,16] for simulating wave propagation are presented in literature. For instance, multiscale finite element methods (MsFEM) [21,22] and the heterogeneous multiscale methods (HMM) [1,2,3] are proposed in order to discretize the wave equation in a coarse grid. Recently, a class of generalized finite element methods for the wave equation [4,23] has been proposed. This type of methods is based on the idea of localized orthogonal decomposition (LOD) [24] and generalize the traditional finite element method to resolve the multiscale problems with a cheaper cost accurately.
In this research, we focus on the recently developed generalized multiscale finite element method (GMsFEM) [7,13]. The GMsFEM is a generalization of the classical MsFEM [14] in the sense that multiple basis functions can be systematically added to each coarse block. The GMsFEM consists of two stages: the offline and online stages. In the offline stage, a set of (local supported) snapshot functions are constructed, which can be used to capture all fine-scale features of the solution essentially. Then, a model reduction is performed by the use of a well-designed local spectral decomposition, and the dominant modes are chosen to be the multiscale basis functions. All these computations are done before the actual simulations of the solution. In the online stage, with a given source term and boundary conditions, the multiscale basis functions obtained in the offline stage are used to approximate the solution. Some previous works are using GMsFEM for the wave equation based on the second-order formulation of wave equation [8,17] and the wave equation in mixed formulation [11].
The objective of this work is to develop for the wave equation in mixed formulation [19] a new computational multiscale method based on the idea of constraint energy minimization (CEM) proposed in [10]. In order to derive an energy-conserving numerical scheme for the wave equation, we consider a pressure-velocity formulation (i.e. mixed formulation). For spatial discretization, we adopt the idea of CEM-GMsFEM presented in [9,10] and propose a multiscale method for heterogeneous wave propagation and construct multiscale spaces for both, the velocity and the pressure variables. In this research, we show the first-order convergence of the method using CEM-GMsFEM combined with the leapfrog scheme. Numerical results are provided to demonstrate the efficiency of the proposed method. The present CEM-GMsFEM setting allows flexibly adding additional basis functions based on spectral properties of the differential operators. This enhances the accuracy of the method in the presence of high contrast in the media. It is shown that if enough basis functions are selected, the convergence of the method can be shown independently of the contrast. Unfortunately, a high number of basis functions directly influences the computational complexity of the method. The direct influence of the contrast on the needed number of basis functions is not known, but numerical results indicate that a moderate number of basis functions, depending logarithmically on the contrast, seems sufficient.
The remainder of the paper is organized as follows. We provide in Section 2 the background knowledge of the problem. Next, we introduce the multiscale method and the discretization in Section 3. In Section 4, we provide the stability estimate of the method and prove the convergence of the proposed method. We present the numerical results in Section 5. Finally, we give some concluding remarks in Section 6.

Preliminaries
Consider the wave equation in mixed formulation over the (bounded) computational domain Ω ⊂ R 2 Here,v andṗ represent the time derivatives of v and p respectively, T > 0 is a given terminal time, ρ ∈ L ∞ (Ω) is the (positive) density of the fluid satisfying 0 < ρ min ≤ ρ, and n is the unit outward normal vector to the boundary ∂Ω. We assume that the permeability field κ : Ω → R is highly oscillatory, satisfying 0 < κ min ≤ κ(x) ≤ κ max for almost every x ∈ Ω with κmax κ min ≫ 1. The source function satisfies f ∈ L 2 (Ω). Here, h v and h p are some given initial conditions. In general, we refer to the solution v as velocity and p as pressure. We denote V 0 := {v ∈ H(div, Ω) : v · n = 0 on ∂Ω} and Q := L 2 (Ω).
Instead of the original PDE formulation, we consider the variational formulation corresponding to (1): find v ∈ V 0 and p ∈ Q such that where (·, ·) denotes the inner product in L 2 (Ω) and (·, ·) ρ denotes the inner product in L 2 (Ω) with weighted function ρ. The bilinear forms a : for all v, w ∈ V 0 and p ∈ Q. We remark that the condition Ω p dx = 0 should be imposed for the pressure variable p ∈ Q in order to make the variational formulation well-posed. Also, the following inf-sup condition should satisfy: for all q ∈ Q with Ω q dx = 0, there exists a constant C 0 > 0 independent to κ such that In this research, we will apply the constraint energy minimizing generalized multiscale finite element method (CEM-GMsFEM) for mixed formulation, which is originally proposed in [10], to approximate the solution of the above mixed problem. First, we introduce fine and coarse grids for the computational domain. Let T H = {K i } N i=1 be a conforming partition of the domain Ω with mesh size H > 0 defined by We refer to this partition as the coarse grid. We denote the total number of coarse elements as N ∈ N + . Subordinate to the coarse grid, we define the fine grid partition T h (with mesh size h ≪ H) by refining each coarse element K ∈ T H into a connected union of finer elements. We assume that the refinement above is performed such that T h is also a conforming partition of the domain Ω. Denote N c as the number of interior coarse grid nodes of T H and we write {x i } Nc i=1 as the interior coarse nodes in the coarse grid T H . The mixed wave problem (1) can be numerically solved on the fine grid T h by the lowest order Raviart-Thomas (RT 0) finite element method. Let (V h , Q h ) be the RT 0 finite element spaces with respect to T h . The approximated variational formulation reads: find We remark that the solution pair (v h , p h ) ∈ V h × Q h is served as a reference solution. In the following sections, we will construct multiscale solution (v ms , p ms ) that gives good ap proximation of (v h , p h ) and derive the corresponding error estimation. For an error bound of the reference solution (v h , p h ), one can apply the technique in [6] where C > 0 is a constant depending on the regularity of the exact solution (v, p).

Methodology
In this section, we outline the framework of CEM-GMsFEM and introduce the construction of the multiscale spaces for approximating the fine-scale solution (v h , p h ). We emphasize that the multiscale basis functions and the corresponding spaces are defined with respect to the coarse grid T H . The multiscale method consists of two steps. First, we construct a multiscale space Q ms for approximating the pressure. Based on the space Q ms , we construct another multiscale space V ms for the velocity. We remark that these basis functions are locally supported in some coarse patches formed by some coarse elements.
Once the multiscale spaces are ready, one can discretize time derivatives in the problem by finite difference and solve the resulting fully discretized problem.

The multiscale method
First, we introduce some notations used in the following sections. Given a subset S ⊂ Ω, where n S is the unit outward normal vector with respect to the boundary ∂S.
The multiscale solution (v ms , p ms ) ∈ V ms × Q ms is obtained by solving the variational formulation We will detail the constructions for the multiscale spaces in next sections.

Construction of pressure basis functions
In this section, we present the construction of the multiscale space Q ms for pressure. For each coarse element K i ∈ T H , consider the following local spectral problem over K i : find for j = 1, · · · , L i , where L i ∈ N + is a local parameter depending on the grids T H and T h .
Here, the bilinear form s i : and {χ j } Nc j=1 is a set of standard multiscale partition of unity. In particular, given an interior coarse grid node x j , the function χ j is defined as the solution to the following system over the coarse neighborhood ω j : where g j is a linear continuous function on all edges of ∂K. Assume that s i (p i j , p i j ) = 1 and we arrange the eigenvalues in ascending order such that 0 ≤ λ 1 corresponding the first J i smallest eigenvalues and we define the multiscale space Q ms for pressure as follows Q ms := span{p i j : i = 1, · · · , N, j = 1, · · · , J i }.

Construction of velocity basis functions
In this section, we present the construction of the multiscale space V ms for velocity. To define the velocity basis, we introduce the operator π : Q h → Q ms as follows Next, we define the bilinear form s : Note that the operator π : Q h → Q ms is the projection of Q h onto the multiscale space Q ms with respect to the inner product s(·, ·). We denote the norm induced by this inner product s(·, ·) as · s .
For a given coarse element K i ∈ T H and ℓ ∈ N + , we define K i,ℓ to be the oversampled region obtained by enlarging ℓ layers from K i . Specifically, we have For simplicity, we denote the oversampled region as K + i . For each eigenfunction p i j ∈ Q ms obtained from (8)-(9), we compute the multiscale basis for velocity ψ i j,ms ∈ V h,0 (K + i ) by solving the following system Then, the multiscale space for velocity is defined as follows V ms := span{ψ i j,ms : i = 1, · · · , N, j = 1, · · · , J i }.

Discretizations
The multiscale spaces V ms and Q ms obtained in Section 3.1 are constructed in the spirit of CEM-GMsFEM. To simplify the notation, we assume that Then, we obtain the following matrix representation from the variational formulation where 0 ∈ R M is the zero vector in R M , ∈ R M are the vectors of coefficients for the approximations v ms and p ms . More precisely, we have For the time discretization, we simply replace the continuous time derivatives by the forward difference in time with a given time step τ > 0. In particular, the velocity term will be approximated at t n = nτ and the pressure term will be approximated at t n+ 1 2 = n + 1 2 τ for n ∈ {0, 1, · · · , N T } with T = N T τ . We remark that the time step τ will be chosen such that N T ∈ N + . It leads to the following fully discretized system where v n = v(t n ), p n+ 1 2 = p(t n+ 1 2 ) and f n := f (t n ), p i M i=1 . We remark that the stability estimate for τ can be obtained by standard techniques and the inverse estimate, see for example [18].

Stability and convergence analysis
In this section, we present the results of stability and convergence analysis for the mixed CEM-GMsFEM established in Section 3. To this aim, we first introduce notation that will be used in this section: We remark that the norms · ρ and · s are equivalent to the standard L 2 -norm · L 2 (Ω) . For p ∈ Q h , we have One may show the equivalence between · s and · L 2 (Ω) provided κ ∈ [κ min , κ max ]. Further, we denote a b if there is a generic constant C > 0 such that a ≤ Cb.

Energy conservation and stability
In this section, we prove the property of energy conservation and the stability of the proposed scheme (6)- (7). In particular, we show the following proposition. Moreover, the following estimate holds Proof. Take v = v ms in (6) and q = p ms in (7). Adding two equations, we obtain a(v ms , v ms ) + (ṗ ms , p ms ) ρ = 1 2 Using the Cauchy-Schwarz inequality, we have It implies that v ms (t, ·) 2 a + p ms (t, ·) 2 ρ ≤ 2 max Using Young's and Jensen's inequalities, we obtain v ms (t, ·) 2 a + p ms (t, ·) 2 ρ ≤ This completes the proof.
Remark. The technique showing (18) from (17) in the proof above will be employed in the convergence analysis below.

Convergence analysis
In this section, we show the convergence result of the semi-discretized scheme. We define (σv, σp) ∈ V ms × Q ms to be the multiscale projection of a given pair ( where Λ := min 1≤i≤N λ i j and {λ i j } are the eigenvalues obtained from (8)- (9).
and the following estimate holds using the result of [10, Theorem 1]. Therefore, we have On the other hand, since b(v − σv, q) = 0 for all q ∈ Q ms , there exists a set of real number Then, by the orthogonality of the eigenfunctions {p i j } and (8)-(9), we have This completes the proof of the lemma.
The main result in this research reads as follows.
Theorem 4.2. Suppose that (v h , p h ) is the solution to the system (4)- (5) and (v ms , p ms ) is the solution to (6)- (7). Then, the following estimate holds where (σv h , σp h ) ∈ V ms × Q ms is the multiscale projection of (v h , p h ) and F := f − ρṗ h .
Rewriting the system above, it implies that Take w = σv h − v ms ∈ V ms and q = σp h − p ms ∈ Q ms in the system above. Then, adding two equations together, we obtain the following equality Here, the properties of multiscale projection (19) are used to simplify the expression above. Hence, we obtain by Cauchy-Schwarz inequality

It implies that
Using the same technique as that of proving Proposition 4.1, we obtain Next, we analyze the terms σ v h −v h 2 a and σ p h −ṗ h 2 ρ . Note that the fine-scale solution (v h , p h ) satisfies the following system where F = f − ρṗ h . Using the properties of the multiscale projection (19), we obtain Then, by the error estimate in [10, Theorem 1], one may obtain the following estimate Similarly, using the same estimate, one may obtain the estimate forσ p h −ṗ h as follows Combining (22), (23), and (24) yields the desire estimate. This completes the proof.
Remark. From the result (20) and the inequality (21), one may easily conclude that v h − v ms a = O(H). Moreover, using the technique in [11,Theorem 5.4] one may show that p h satisfies

Numerical experiments
In this section, we perform some numerical experiments using the proposed multiscale method to solve the wave equation in mixed formulation. Let the computational domain and the density be Ω = (0, 1) 2 , T = 0.4 and ρ ≡ 1. In the simulation, we use (uniform) rectangular mesh to perform the spatial discretization and set the fine mesh size to be h = √ 2/400. We will specify the coarse mesh size H > 0 in the examples below. For the time discretization, set the time step to be τ = 10 −4 . In both examples, let the initial conditions be h v = 0 and h p = 0. Recall that ℓ ∈ N is the number of oversampling layers used to perform the constraint energy minimization. We set J := J i as the number of basis functions to from the auxiliary space Q ms .
We use the following quantities of error to measure the performance of the proposed scheme where (v h , p h ) is the fine-scale solution computed on T h and (v ms , p ms ) is the multiscale solution obtained by solving (6)-(7).
Example 5.1 (Heterogeneous model). In this example, we consider the case with highcontrast permeability field κ over the computational domain. See Figure 1(a) for an illustration of this permeability field. The source function f : R 2 → R is set to be Clearly, it holds that f ∈ L 2 (Ω) and it satisfies Ω f dx = 0. The solution profiles of the  We explore the efficiency of the proposed method by showing the errors in velocity and pressure. Tables 1 and 2 show the quantities e vel and e pre in different time levels with different H while fixing J = 5 and ℓ = 3. One can see that the error (either in pressure or velocity) decreases as the coarse mesh size decreases.      Next, we fix the coarse mesh size H = √ 2/20 and the number of oversampling layers ℓ = 4. We adjust the number of basis functions J to see how this factor affects the errors in velocity and pressure. Tables 3 and 4 show the errors e vel and e pre . Under this setting, one may observe that the errors e pre and e vel are reduced as the number of basis functions J increases. We remark that once the number of basis functions reaches a certain level, the decay of error becomes slower due to the fact that the decay of eigenvalues obtained in (8)    layers. Tables 5 and 6 show the corresponding results. One may observe that the accuracy of the solution is improved if more layers are included in the simulation. Once the number of layers ℓ exceeds a certain level, the decay of error stagnates. We remark that based on the theoretical findings in [10], the number of layers should depend on the logarithm of the maximum value of the contrast.    Figure 1(b).
The source function f is chosen as the first derivative of the Gaussian wavelet with central for x ∈ Ω, t ∈ (0, T ] and c = (0.5, 0.5) T . Note that δ > 0 measures the size of the support of the source. Here, we denote the two-dimensional Euclidean distance as |x − y| for any x, y ∈ R 2 .   Similar as before, we report the results of errors in both velocity and pressure using the proposed multiscale method with varying coarse mesh size H. Tables 7 and 8 record the errors with fixed J = 5, ℓ = 3, f 0 = 20 and δ = 0.02. We remark that δ = 4 √ 2h in this case. One can observe that the errors in pressure and velocity reduce as the coarse mesh size decreases.
Furthermore, we calculate the multiscale solution by the proposed method with different f 0 and δ when H = √ 2/20, J = 5 and ℓ = 3. Tables 9 and 10 present the corresponding numerical results. We remark that the For high frequency cases with larger f 0 , one may use a smaller H to enhance the accuracy of the multiscale approximation.

Conclusion
In this work, we have proposed and analyzed the constraint energy minimizing generalized multiscale finite element method for solving the wave equation in mixed formulation. The multiscale basis functions for pressure are obtained by solving a class of well-designed local spectral problems. Based on the concept of constraint energy minimization, we construct the multiscale basis functions for velocity satisfying the property of least energy. The method is shown to have first-order convergence with respect to the coarse mesh size. Numerical results are provided to illustrate the efficiency of the proposed method.