TY - JOUR
AU1 - Meng,, Weijuan
AU2 - Fu,, Li-Yun
AB - Abstract The finite element method is a very important tool for modeling seismic wave propagation in complex media, but it usually consumes a large amount of memory which significantly decreases computational efficiency when solving large-scale seismic problems. Here, a modified finite element method (MFEM) is proposed to improve efficiency. Triangular elements are employed to mesh the topography and the discontinuous interface more flexibly. In the two-dimensional case, the Jacobian matrix is obtained by using three controlling points instead of all nodes in each element with MFEM, which separates the Jacobian matrix from the stiffness matrix. The kernel matrices of the stiffness matrix rather than the global matrix are stored, and memory requirements are thus reduced significantly. Meanwhile, the element-by-element scheme is adopted to spare large sparse matrices and make the program easily parallelized. A second-order perfectly matched layer (PML) is also implemented to eliminate artificial reflections. Finally, the accuracy and efficiency of our algorithm are validated by numerical tests. numerical simulation, finite element method, storage scheme, element-by-element, PML absorbing boundary condition 1. Introduction Investigation of the wavefield from the earth’s interior as well as the interpretation of seismic data networks have been major problems. Therefore, the numerical simulation of seismic wave propagation is very important in the seismology (Kelly et al1976, Komatitsch and Tromp 2002). Ongoing computer hardware improvements in recent decades have provided the possibility of extracting more information from real seismic data. Widely used methods for solving seismic wave equations include the finite difference method (FDM) (Alford et al1974, Virieux 1984, 1986, Bohlen and Saenger 2006), the pseudo spectral method (PSM) (Fornberg 1988, Liu et al2013), the finite element method (FEM) (Padovani et al1994, Liu et al2014a), the spectral element method (SEM) (Komatitsch and Vilotte 1998, Komatitsch and Tromp 2002, 2003, Komatitsch et al2010) and other optimized methods (Yang et al2009). The FDM is most widely used due to its high efficiency and easy implementation. However, the conventional FDM has some disadvantages in practical modeling. Firstly, severe numerical dispersion occurs when a coarse grid is used or a high velocity contrast exists. Secondly, rugged surfaces are hard to subdivide directly and a free surface boundary condition cannot be satisfied naturally and needs to be dealt with by extra equations and other special schemes (Kelly et al1976, Moczo et al2000, Li et al2010, 2011, 2012, Zhang and Shen 2010, Zhang et al2012, Wang et al2014). The PSM is a global method with good precision and little numerical dispersion. However, it may yield numerical instability in a strongly heterogeneous medium and is still incapable of handling complex geological geometries or free boundary conditions (Fornberg 1988, Tessmer and Kosloff 1994, Liu et al2013). The FEM overcomes some shortcomings of previous methods. It is applicable to various geometrical models and naturally satisfies the free boundary condition. However, the large memory requirements decrease its computational efficiency (Marfurt 1984, Padovani et al1994, Zhang et al2002, Liu et al2014a). The SEM, known as the high-order FEM, also needs a large amount of computer memory. One benefit of SEM with Legendre polynomials is that the interpolation points are consistent with the integral points, which leads to a diagonal mass matrix. Another is that a spectral-like convergence rate can be achieved because the orthogonal polynomials are used in the space domain. Therefore, it has high computational accuracy and efficiency in modeling seismic wave (Komatitsch and Vilotte 1998, Chaljub et al2007, Komatitsch et al2010). However, the rectangular elements (two-dimensional, 2D) and hexahedral elements (three-dimensional, 3D) used in conventional SEM can hardly subdivide the complex model (Liu et al2014). In the classical FEM either the consistent mass matrix or the stiffness matrix takes up most computer memory and decreases computational efficiency dramatically. The original mass matrix, also known as the consistent mass matrix, is non-diagonal and consumes a large amount of computer memory (Marfurt 1984, Padovani et al1994). Based on the law of mass conservation, the whole mass of each element is dispersed along its nodes. Then, the consistent mass matrix is transformed into a lumped mass matrix, which is a diagonal matrix (Richter 1994, Liu et al2014b). Such special treatment of the mass matrix avoids the need to solve a large-scale sparse system of equations and substantially increases the computational efficiency of the FEM. Each discrete point in the computational domain of the FEM is only connected with the points in the surrounding elements. Thus, most zero elements are in each row of the stiffness matrix and the non-zero ones tend to be distributed around the main diagonal of the stiffness matrix (Jiang and Pang 1979). Constant bandwidth storage is a simple and common way to store the stiffness matrix (Cuthill and McKee 1969, Gibbs et al1976). In this scheme, only the elements inside the semi-bandwidth of the stiffness matrix are stored and the zero elements outside can be neglected. The drawback of this storage mode is that a large number of zero elements are stored inside the semi-bandwidth. In order to reduce storage requirements, some data-compression storage techniques, such as the compressed sparse row (CSR) method (Liu et al2013), are introduced. In the CSR method, only non-zero elements in each row of the stiffness matrix are stored in three one-dimensional arrays. Therefore, computer memory and computational cost can both be reduced. All these previous methods cannot avoid storing the large sparse stiffness matrix. Thus, the storage costs of these methods may still be unaffordable in solving large-scale problems. Liu et al (2014a) introduced the kernel matrices storage (KMS) method in which only the kernel matrices of the global matrix are stored and the assembly of a large global stiffness matrix can be avoided. The matrix–vector products are calculated element-by-element, similar to the discontinuous Galerkin methods (Hesthaven and Warburton 2008, He et al2015). At the cost of a small increase in float-point computation time, the storage memory for the KMS method can be significantly reduced. However, KMS is only applicable to consistent elements and is less effective in complex models with irregular interfaces, and thus requires further improvement. This paper is organized as follows. In the next section, the improved version of KMS method, called the modified finite element method (MFEM), is proposed. The lumped mass matrix, the second-order perfectly matched layer (PML) absorbing boundary condition and updating of the wavefield are also presented. Then the accuracy and efficiency of the MFEM are investigated in a simple homogeneous medium by comparison with different storage schemes and other classical methods such as conventional FEM and SEM. In addition, the effectiveness of different-order MFEMs are also benchmarked in a two-layered model with inclined interfaces and a rugged surface model. Finally, the computational accuracy and efficiency of the MFEM are discussed. 2. Theoretical method 2.1. Elastic wave equation The elastic wave equation is used to describe the propagation of a seismic wave through the solid earth. Considering an isotropic elastic medium, the 2D elastic wave equation is given by (Ma et al 2011, Liu et al2015b) ρ∂2u∂t2=∇((λ+μ)∇⋅u)+∇⋅(μ∇u)+f,1 where u and f are the displacement vector and source term, respectively, ρ is mass density and λ and μ are Lamé constants. Multiplying the above equation by a test function, integrating over the whole domain Ω and applying the method of weighted residuals, the weak form of equation (1) yields ∫Ωρ∂2u∂t2⋅vdx=∫Ω(λ+μ)(∇⋅u)(∇⋅v)dx+∫Ωμ(∇u):(∇v)dx+∫Ωf⋅vdx,2 where v is an arbitrary time-independent vector and x is the Cartesian coordinate. We subdivide the domain Ω into N non-overlapping elements, and equation (2) can be written as MÜ=KU+FM=∑eN∫Ωeρ[φ]T[φ]dΩK=∑eNK11K12K21K22ΩeK11=-∫Ωe(λ+2μ)[φx]T[φx]+μ[φz]T[φz]dΩK12=-∫Ωeλ[φx]T[φz]+μ[φz]T[φx]dΩK21=-∫Ωeμ[φx]T[φz]+λ[φz]T[φx]dΩK22=-∫Ωeμ[φx]T[φx]+(λ+2μ)[φz]T[φz]dΩ3 where U is the displacement vector, M is the consistent mass matrix, K is the assembled stiffness matrix and F is the source vector. For regular quadrilateral elements, each stiffness kernel matrix is identical in all the same shape elements, and only four kernels, ∫Ωe[φx]T[φx]dΩ, ∫Ωe[φz]T[φz]dΩ, ∫Ωe[φx]T[φz]dΩ and ∫Ωe[φz]T[φx]dΩ, of an element stiffness matrix instead of a huge global matrix are stored. A detailed description can be found in Liu et al (2014a). However, in order to simulate seismic wave propagation in complex structures we need to investigate this problem with triangular subdivision further. Different triangular elements have different Jacobian matrices that connect the physical coordinates with the local coordinates. Therefore, we only need to calculate Jacobian matrices over all elements and store them in advance. The physical coordinates (x,z) and isoparametric coordinates (ξ,η) of an arbitrary element, e, are shown in figure 1. Figure 1. Open in new tabDownload slide The physical coordinates (x,z) and isoparametric coordinates (ξ,η) of a triangular element e. Figure 1. Open in new tabDownload slide The physical coordinates (x,z) and isoparametric coordinates (ξ,η) of a triangular element e. Obviously, the three vertices, X1, X2 and X3, in area coordinates are A(1, 0), B(0, 1) and C(0, 0). According to the relation between these two coordinate systems, we can easily obtain the following expression x=(x1-x3)ξ+(x2-x3)η+x3z=(z1-z3)ξ+(z2-z3)η+z3.4 Here the relation between the coordinates (x,z) and (ξ,η) is linear. Therefore, the Jacobian matrix J, which maps the partial derivative in terms of physical coordinates (x,z) to the partial derivative in terms of area coordinates (ξ,η), can be obtained from linear interpolation ∂∂ξ∂∂η=J∂∂x∂∂z=∂x∂ξ∂z∂ξ∂x∂η∂z∂η∂∂x∂∂z.5 The elements of J can be obtained as follows ∂x∂ξ=∑i=13xiφξi(ξ,η)=(x1-x3)∂z∂ξ=∑i=13ziφξi(ξ,η)=(z1-z3)∂x∂η=∑i=13xiφηi(ξ,η)=(x2-x3)∂z∂η=∑i=13ziφηi(ξ,η)=(z2-z3).6 Since the Jacobian matrix is fixed for a particular element, it can be calculated using only three controlling points instead of all points in the element. Different triangular elements have different Jacobian matrices. Thus, the stored kernel matrices should be free of the impact of the Jacobian matrix. According to the transformation principle of two coordinate systems, we have φx=∂ξ∂xφξ+∂η∂xφηφz=∂ξ∂zφξ+∂η∂zφη7 and ∫Ωe[φx]T[φx]dΩ=∂ξ∂x∂ξ∂x∫Ωe[φξ]T[φξ]dΩ+∂ξ∂x∂η∂x∫Ωe[φξ]T[φη]dΩ+∂η∂x∂ξ∂x∫Ωe[φη]T[φξ]dΩ+∂η∂x∂η∂x∫Ωe[φη]T[φη]dΩ∫Ωe[φx]T[φz]dΩ=∂ξ∂x∂ξ∂z∫Ωe[φξ]T[φξ]dΩ+∂ξ∂x∂η∂z∫Ωe[φξ]T[φη]dΩ+∂η∂x∂ξ∂z∫Ωe[φη]T[φξ]dΩ+∂η∂x∂η∂z∫Ωe[φη]T[φη]dΩ∫Ωe[φz]T[φx]dΩ=∂ξ∂z∂ξ∂x∫Ωe[φξ]T[φξ]dΩ+∂ξ∂z∂η∂x∫Ωe[φξ]T[φη]dΩ+∂η∂z∂ξ∂x∫Ωe[φη]T[φξ]dΩ+∂η∂z∂η∂x∫Ωe[φη]T[φη]dΩ∫Ωe[φz]T[φz]dΩ=∂ξ∂z∂ξ∂z∫Ωe[φξ]T[φξ]dΩ+∂ξ∂z∂η∂z∫Ωe[φξ]T[φη]dΩ+∂η∂z∂ξ∂z∫Ωe[φη]T[φξ]dΩ+∂η∂z∂η∂z∫Ωe[φη]T[φη]dΩ.8 Here, the kernel matrices are ∫Ωe[φξ]T[φξ]dΩ, ∫Ωe[φξ]T[φη]dΩ, ∫Ωe[φη]T[φξ]dΩ and ∫Ωe[φη]T[φη]dΩ. From the above discussion, it can be inferred that the global stiffness matrix can be easily obtained if the Jacobian matrix of triangular elements as well as the kernel matrices are calculated in advance. 2.2. Lumped mass matrix The consistent mass matrix constructed from the classical finite element equations is non-diagonal and requires the solution of a large sparse linear system of equations. However, based on the principle of mass conservation, the non-diagonal consistent mass matrices can be lumped into a diagonal mass matrix with the lumped mass technique (Archer and Whalen 2005, Wu 2006). The formula for diagonalizing the mass matrix M in equation (3) is given by (Liu et al2014b) Miid=Mii∑i∑jMij/∑jMjj.9 Because of the diagonal property of the lumped mass matrix, the wavefields can be updated by explicit recursion and the efficiency of the FEM can be improved significantly. However, when the mesh size is much less than the seismic wavelength and the accuracy of numerical integration is adequate, the numerical precision of the FEM with the lumped mass matrix is of the same order as that with the consistent mass matrix. Here, we do not discuss the numerical precision specifically and just test the accuracy of our method by comparing it with analytical solutions or other methods. More details can be found in the papers of Wu (2006) or Liu et al (2013). 2.3. PML absorbing boundary condition In order to eliminate artificial reflections, the PML for second-order elastic wave equations is adopted (Komatitsch and Tromp 2003, Liu et al2012). From equation (1), the source-free elastic wave equation can be rewritten in the following form: ρ∂2ux∂t2=∂∂x(λ+2μ)∂ux∂x+λ∂uz∂z+∂∂zμ∂uz∂x+∂ux∂zρ∂2uz∂t2=∂∂xμ∂uz∂x+∂ux∂z+∂∂zλ∂ux∂x+(λ+2μ)∂uz∂z.10 In complex stretched coordinates (Collino and Tsogka 2001), equation (10) can be rewritten as ρ(iω)2Ux=iωiω+dx(x)∂∂x(λ+2μ)iωiω+dx(x)×∂Ux∂x+λiωiω+dz(z)∂Uz∂z+iωiω+dz(z)∂∂zμiωiω+dx(x)∂Uz∂x+iωiω+dz(z)∂Ux∂zρ(iω)2Uz=iωiω+dx(x)∂∂xμiωiω+dx(x)∂Uz∂x+iωiω+dz(z)∂Ux∂z+iωiω+dz(z)∂∂zλiωiω+dx(x)∂Ux∂x+(λ+2μ)iωiω+dz(z)∂Uz∂z11 where i=-1 is the imaginary unit, ω is the angular frequency and Ux and Uz are the Fourier transforms of ux and uz with respect to time, respectively. dx(x) and dz(z) are the damping coefficients in the PML region, which are described in detail in the appendix. After splitting each displacement component into three terms and transforming the terms back into the time domain, we get ρ[∂t+dx(x)]2ux1=∂∂x(λ+2μ)∂ux∂x+ρPxxρ[∂t+dz(z)]2ux2=∂∂zμ∂ux∂z+ρPxzρ[∂t+dx(x)][∂t+dz(z)]ux3=∂∂xλ∂uz∂z+∂∂zμ∂uz∂xρ[∂t+dx(x)]Pxx=-(λ+2μ)dx′(x)∂ux∂xρ[∂t+dz(z)]Pxz=-μdz′(z)∂ux∂z12 ρ[∂t+dz(z)]2uz1=∂∂z(λ+2μ)∂uz∂z+ρPzzρ[∂t+dx(x)]2uz2=∂∂xμ∂uz∂x+ρPzxρ[∂t+dx(x)][∂t+dz(z)]uz3=∂∂xμ∂uz∂z+∂∂zλ∂uz∂xρ[∂t+dz(z)]Pzz=-(λ+2μ)dz′(z)∂uz∂zρ[∂t+dx(x)]Pzx=-μdx′(x)∂uz∂x13 where ux=ux1+ux2+ux3 and uz=uz1+uz2+uz3. Pxx, Pxz, Pzx and Pzz are four auxiliary variables. Multiplying equations (12) and (13) with an arbitrary test function, v=(vx,vz)T, and conducting partial integration, we can obtain the weak forms: ∫Ω2ρ[∂t+dx(x)]2ux1vxdΩ=-∫Ω2(λ+2μ)×∂ux∂x∂vx∂xdΩ+∫Ω2ρPxxvxdΩ∫Ω2ρ[∂t+dz(z)]2ux2vxdΩ=-∫Ω2μ∂ux∂z∂vx∂zdΩ+∫Ω2ρPxzvxdΩ∫Ω2ρ[∂t+dx(x)][∂t+dz(z)]ux3vxdΩ=-∫Ω2λ∂uz∂z∂vx∂xdΩ-∫Ω2μ∂uz∂x∂vx∂zdΩ∫Ω2ρ[∂t+dx(x)]PxxvxdΩ=-∫Ω2(λ+2μ)dx′(x)∂ux∂xvxdΩ∫Ω2ρ[∂t+dz(z)]PxzvxdΩ=-∫Ω2μdz′(z)∂ux∂zvxdΩ,14 ∫Ω2ρ[∂t+dz(z)]2uz1vzdΩ=-∫Ω2(λ+2μ)×∂uz∂z∂vz∂zdΩ+∫Ω2ρPzzvzdΩ∫Ω2ρ[∂t+dx(x)]2uz2vzdΩ=-∫Ω2μ∂uz∂x∂vz∂xdΩ+∫Ω2ρPzxvzdΩ∫Ω2ρ[∂t+dx(x)][∂t+dz(z)]uz3vzdΩ=-∫Ω2μ∂uz∂z∂vz∂xdΩ-∫Ω2λ∂uz∂x∂vz∂zdΩ∫Ω2ρ[∂t+dz(z)]PzzvzdΩ=-∫Ω2(λ+2μ)dz′(z)∂uz∂zvzdΩ∫Ω2ρ[∂t+dx(x)]PzxvzdΩ=-∫Ω2μdx′(x)∂uz∂xvzdΩ,15 where Ω2 denotes the PML domain. After discretizing equations (14) and (15), we can obtain the finite element equations as follows: M∂t2Ux1+2Cx∂tUx1=MPxx-Kx1Ux-CxxUx1M∂t2Ux2+2Cz∂tUx2=MPxz-Kx2Ux-CzzUx2M∂t2Ux3+(Cx+Cz)∂tUx3=-Kx3Uz-CxzUx3M∂tPxx+CxPxx=-Kx4UxM∂tPxz+CzPxz=-Kx5Ux,16 M∂t2Uz1+2Cz∂tUz1=MPzx-Kz1Uz-CxxUz1M∂t2Uz2+2Cx∂tUz2=MPzz-Kz2Uz-CzzUz2M∂t2Uz3+(Cx+Cz)∂tUz3=-Kz3Ux-CxzUz3M∂tPzx+CxPzx=-Kz4UzM∂tPzz+CzPzz=-Kz5Uz,17 where Ux=Ux1+Ux2+Ux3 and Uz=Uz1+Uz2+Uz3. The coefficient matrices are expressed as: M=∑Ω2∫Ωeρ[φ]T[φ]dΩ,Cx=∑Ω2∫Ωeρdx(x)[φ]T[φ]dΩCz=∑Ω2∫Ωeρdz(z)[φ]T[φ]dΩ,Cxx=∑Ω2∫Ωeρdx2(x)[φ]T[φ]dΩCzz=∑Ω2∫Ωeρdz2(x)[φ]T[φ]dΩ,Cxz=∑Ω2∫Ωeρdx(x)dz(z)[φ]T[φ]dΩKx1=∑Ω2∫Ωe(λ+2μ)[φx]T[φx]dΩ,Kx2=∑Ω2∫Ωeμ[φz]T[φz]dΩKx3=∑Ω2∫Ωe{λ[φx]T[φz]+μ[φz]T[φx]}dΩ,Kx4=∑Ω2∫Ωe(λ+2μ)[dx(x)][φ]T[φx]T[φx]dΩKx5=∑Ω2∫Ωeμ[dz(z)][φ]T[φz]T[φz]dΩ,Kz1=∑Ω2∫Ωe(λ+2μ)[φz]T[φz]dΩKz2=∑Ω2∫Ωeμ[φx]T[φx]dΩ,Kz3=∑Ω2∫Ωe{μ[φx]T[φz]+λ[φz]T[φx]}dΩKz4=∑Ω2∫Ωeμ[dx(x)][φ]T[φx]T[φx]dΩ,Kz5=∑Ω2∫Ωe(λ+2μ)[dz(z)][φ]T[φz]T[φz]dΩ.18 2.4. Updating the wavefield There are many methods for temporal discretization, such as the Lax–Wendroff method (Chen 2009), the Runge–Kutta method (Yang et al2009), the Newmark method (Komatitsch and Tromp 2002, Liu et al2014) and symplectic schemes (Ma et al 2011, Liu et al2015a). Among these methods, the second-order Lax–Wendroff method, also known as the central difference method, is employed due to its simplicity and high efficiency: U(t+Δt)=(Δt)2M-1(KU(t)+F)+2U(t)-U(t-Δt).19 3. Numerical experiments 3.1. A simple homogeneous model First, we use a homogeneous medium to investigate the validity and the accuracy of the MFEM. The size of the model is 1000 m × 1000 m, and the thickness of the PML absorbing boundary is 100 m. The velocities of P- and S-waves are 2000 m s-1 and 1250 m s-1, respectively. The density of the material is 2000 kg m-3. A vertical force is located at (500 m, -500 m). The source time function is a Ricker wavelet with a center frequency of 25 Hz. The time increment is 0.2 ms. The geometry of the physical model is shown in figure 2(a). In order to compare our method with the same order SEM, we discretize the model into structured triangular elements with a size of 10 m. The schematic diagram for subdivision of the computational domain is shown in figure 2(b). Figure 2. Open in new tabDownload slide (a) The geometry of an elastic homogeneous model. The source is located at (500 m, -500 m). The two receivers R1 and R2 are placed at (500 m, 0 m) and (800 m, 0 m), respectively. (b) Schematic diagram of the structured triangular subdivision. Figure 2. Open in new tabDownload slide (a) The geometry of an elastic homogeneous model. The source is located at (500 m, -500 m). The two receivers R1 and R2 are placed at (500 m, 0 m) and (800 m, 0 m), respectively. (b) Schematic diagram of the structured triangular subdivision. In addition, we present the memory consumptions of the stiffness matrix in different-order FEMs and various storage schemes in table 1. M1, denoting a FEM with all elements of the global stiffness matrix stored, requires a huge memory and cannot be implemented on any computer. For constant bandwidth storage, the node number has a great effect on the bandwidth and the memory requirement. Thus, M2 is presented as the lower limit to the constant bandwidth storage scheme. M3 denotes FEM with CSR and M4 denotes MFEM with KMS. Taking second-order methods as an example, the number of elements is n1 = 20 000 and the number of nodes including controlling and interpolating nodes is n2 = 40 401. For M1, the number of storage elements is n2 × n2 × 4 = 6528 963 204 and the memory cost is 6528 963 204 × 4 ≈ 249 06.0181 MB in single precision storage. For M2, the bandwidth of the stiffness matrix is nb = 19. The number of storage elements is nb × n2 × 4 = 3070 476 and the memory cost is 3070 476 × 4 ≈ 11.7129 MB. For M3, the numbers of the rows that include six, nine and nineteen non-zero elements are 402, 29 802, 396 and 9801. Thus, the number of storage elements, including non-zero elements over all rows, is 1846 404, and the memory cost for storing these non-zero ones by using the three arrays is 1846 404 × (4 + 1) + (n2 + 1) × 4 ≈ 8.9585 MB. For M4, the number of nodes in each element is n = 6. The number of storage elements is (n1 × 4) + (n × n × 4) =80 144 and the memory cost is 80 144 × 4 ≈ 0.3057 MB. Compared with the memory requirements of constant bandwidth storage and CSR, M4 can save a significant amount of memory. Table 1. Memory requirements of the FEM with different storage schemes. Methods . M1 . M2 . M3 . M4 . First-order Number of storage elements 416 241 604 285 628 282 404 80 036 Memory (MB) 1587.8357 1.0896 1.3855 0.3053 Second-order Number of storage elements 6528 963 204 3070 476 1846 404 80 144 Memory (MB) 24 906.0181 11.7129 8.9585 0.3057 Third-order Number of storage elements 32 834 164 804 13 408 948 6132 004 80 400 Memory (MB) 125 252.3987 51.1511 29.5853 0.3067 Fourth-order Number of storage elements 103 427 846 404 39 235 444 15 059 204 80 900 Memory (MB) 394 545.9229 149.6713 72.4213 0.3086 Methods . M1 . M2 . M3 . M4 . First-order Number of storage elements 416 241 604 285 628 282 404 80 036 Memory (MB) 1587.8357 1.0896 1.3855 0.3053 Second-order Number of storage elements 6528 963 204 3070 476 1846 404 80 144 Memory (MB) 24 906.0181 11.7129 8.9585 0.3057 Third-order Number of storage elements 32 834 164 804 13 408 948 6132 004 80 400 Memory (MB) 125 252.3987 51.1511 29.5853 0.3067 Fourth-order Number of storage elements 103 427 846 404 39 235 444 15 059 204 80 900 Memory (MB) 394 545.9229 149.6713 72.4213 0.3086 Open in new tab Table 1. Memory requirements of the FEM with different storage schemes. Methods . M1 . M2 . M3 . M4 . First-order Number of storage elements 416 241 604 285 628 282 404 80 036 Memory (MB) 1587.8357 1.0896 1.3855 0.3053 Second-order Number of storage elements 6528 963 204 3070 476 1846 404 80 144 Memory (MB) 24 906.0181 11.7129 8.9585 0.3057 Third-order Number of storage elements 32 834 164 804 13 408 948 6132 004 80 400 Memory (MB) 125 252.3987 51.1511 29.5853 0.3067 Fourth-order Number of storage elements 103 427 846 404 39 235 444 15 059 204 80 900 Memory (MB) 394 545.9229 149.6713 72.4213 0.3086 Methods . M1 . M2 . M3 . M4 . First-order Number of storage elements 416 241 604 285 628 282 404 80 036 Memory (MB) 1587.8357 1.0896 1.3855 0.3053 Second-order Number of storage elements 6528 963 204 3070 476 1846 404 80 144 Memory (MB) 24 906.0181 11.7129 8.9585 0.3057 Third-order Number of storage elements 32 834 164 804 13 408 948 6132 004 80 400 Memory (MB) 125 252.3987 51.1511 29.5853 0.3067 Fourth-order Number of storage elements 103 427 846 404 39 235 444 15 059 204 80 900 Memory (MB) 394 545.9229 149.6713 72.4213 0.3086 Open in new tab When the order increases from one to four, the memory cost of the stiffness matrix will increase by 15.69, 78.88 and 248.48 times for M1, by 10.75, 46.95 and 137.37 times for M2, 6.47, 21.35, by 52.27 times for M3 and by about 1 time for M4. The memory cost of the constant bandwidth storage method falls between M1 and M2. Therefore, for the constant bandwidth storage method and M3, the computational points and memory requirements of the stiffness matrix will increase rapidly with increasing numerical order. For M4, only the kernel matrices of elements include more nodes and thus the increase in memory requirements will be trivial. Simulation of seismic wave propagation is performed on the same computer with three numerical methods including the FEM, MFEM, and SEM, which are third-order in space. The computational domain is subdivided into quadrilateral elements in the SEM with a side length of 10 m and triangular elements in the FEM and MFEM with a right side length of 10 m. Figure 3 shows seismic waveforms and errors recorded in the two receivers. The black lines indicate the analytical solutions generated by the Cagniard–De Hoop method (De Hoop 1960), and the other lines indicate the results from FEM, MFEM and SEM. As shown in figure 3, the amplitude and arrival time of the P- and S-waves are accurate, which means that the results calculated using the MFEM and the other methods are consistent with the analytical solutions in the homogeneous model. The maximum norm of errors in FEM, MFEM and SEM are 0.0504, 0.0513 and 0.0207, respectively, at receiver R1 (figure 3(b)). At receiver R2 (figure 3(d)), the errors are 0.0473, 0.0458 and 0.0195, respectively. The results of FEM and SEM with triangular elements are less accurate than those with quadrilateral elements (Liu et al2014). Therefore, the above error analyses are reasonable. Figure 3. Open in new tabDownload slide Seismic records and errors recorded at two receivers: (a), (c) normalized displacement in the z-direction at receivers R1 and R2, respectively; (b), (d) the errors between numerical and analytical solutions. The black lines indicate the analytical solutions calculated with the Cagniard–De Hoop method. The green, red and blue lines indicate the results calculated using FEM, MFEM and SEM, respectively. Figure 3. Open in new tabDownload slide Seismic records and errors recorded at two receivers: (a), (c) normalized displacement in the z-direction at receivers R1 and R2, respectively; (b), (d) the errors between numerical and analytical solutions. The black lines indicate the analytical solutions calculated with the Cagniard–De Hoop method. The green, red and blue lines indicate the results calculated using FEM, MFEM and SEM, respectively. As shown in figure 4(a), the performance and numerical stability of the PML are tested at long times period by measuring the decay of kinetic, potential and total energy versus time. The seismic records of R2 are presented in figures 4(b) and (c). The total energy is defined as (Shi et al2012) E=12∫Ω(ρ∥υ∥2+σ:ε)dΩ.20 We present the evolution of energy over 5 × 104 time steps (10 s). As the source loads over time, the kinetic, potential and total energy significantly increase and remain unchanged from 0.016 to 0.27 s. The energy of the P- and S-waves decreases quickly at about 0.28 s and 0.43 s, at which the waves propagate gradually away from the computational domain. The energy decreases approximately to zero and the seismic records of R2 are close to zero after about 0.64 s, which shows the good performance and stability of the PML without apparent artificial reflection energy at long times. Figure 4. Open in new tabDownload slide (a) The decay of the energies in the inner domain (not including PML zones). The blue, green and red lines represent the kinetic, potential and total energy, respectively. (b), (c) The horizontal and vertical displacements, respectively, recorded by R2. Figure 4. Open in new tabDownload slide (a) The decay of the energies in the inner domain (not including PML zones). The blue, green and red lines represent the kinetic, potential and total energy, respectively. (b), (c) The horizontal and vertical displacements, respectively, recorded by R2. Table 2 shows the computational efficiency of the third- and fourth-order FEM, MFEM and SEM in the homogeneous model. The stiffness matrices are stored with the CSR scheme in FEM and SEM but with the KMS scheme in MFEM. Despite including the computation of PML and other additional costs, MFEM has the smallest memory requirement. In the third-order case, the memory for MFEM is about 0.46 times that for FEM and 0.34 times that for SEM. In the fourth-order case, the memory for MFEM is about 0.41 times that for FEM and 0.28 times that for SEM. In the MFEM, the Jacobian matrices of elements and the kernel matrices are stored ahead and then the equations are resolved element-by-element. Therefore, the MFEM needs extra CPU time (about 10% more than for FEM). The CPU time for MFEM is still less than that for high-order SEM. These conclusions are consistent with the results of Liu et al (2014a), in which only elements of the same shape, namely regular quadrilateral elements, are investigated. Table 2. Computational efficiency of FEM, MFEM and SEM in the first model. Methods . FEM . MFEM . SEM . Third-order Memory (KB) 131 636 61 124 179 564 CPU time (s) 405.51 436.08 457.81 Fourth-order Memory (KB) 304 092 124 948 439 076 CPU time (s) 820.16 926.77 998.99 Methods . FEM . MFEM . SEM . Third-order Memory (KB) 131 636 61 124 179 564 CPU time (s) 405.51 436.08 457.81 Fourth-order Memory (KB) 304 092 124 948 439 076 CPU time (s) 820.16 926.77 998.99 Open in new tab Table 2. Computational efficiency of FEM, MFEM and SEM in the first model. Methods . FEM . MFEM . SEM . Third-order Memory (KB) 131 636 61 124 179 564 CPU time (s) 405.51 436.08 457.81 Fourth-order Memory (KB) 304 092 124 948 439 076 CPU time (s) 820.16 926.77 998.99 Methods . FEM . MFEM . SEM . Third-order Memory (KB) 131 636 61 124 179 564 CPU time (s) 405.51 436.08 457.81 Fourth-order Memory (KB) 304 092 124 948 439 076 CPU time (s) 820.16 926.77 998.99 Open in new tab 3.2. A two-layered model A two-layered model with an inclined interface in the elastic half-space is considered here. The geometry of the model is shown in figure 5. The model is 4000 m wide and 3000 m deep. The inclined interface has a dip of approximately 10°. The upper layer is characterized by a P-wave velocity of 2000 m s-1, an S-wave velocity of 1300 m s-1 and a density of 1000 kg m-3. The lower layer is characterized by a P-wave velocity of 2800 m s-1, an S-wave velocity of 1473 m s-1 and a density of 1500 kg m-3. The source located at (1000 m, -30 m) is a Ricker wavelet with a central frequency of 15 Hz. Figure 5. Open in new tabDownload slide The geometry of the two-layered model. The source is located at (1000 m, -30 m). The two receivers are placed at (1600 m, 0 m) and (1500 m, -1200 m), respectively. Figure 5. Open in new tabDownload slide The geometry of the two-layered model. The source is located at (1000 m, -30 m). The two receivers are placed at (1600 m, 0 m) and (1500 m, -1200 m), respectively. Snapshots of horizontal and vertical displacements at t = 1.0 s are shown in figure 6. Both MFEM and SEM have a good performance in modeling wave propagation in the two-layered model. The snapshots show that all types of waves can be clearly identified and no obvious numerical dispersion occurs. The strong Rayleigh wave indicated by arrow a is observed propagating along the free surface and decays rapidly with depth. The converted S-waves marked by arrow b near the free surface are also clear in the snapshots. The converted waves indicated by arrow c are caused by reflection and transmission of the inclined interface. The PML absorbing boundary attenuates the waves propagating rapidly across the left boundary and thus no artificial reflection occurs. Figure 6. Open in new tabDownload slide Snapshots of wave propagation in a two-layered model with a tilted interface at t = 1.0 s. (a), (b) The horizontal and vertical components of the displacement calculated with MFEM. (c), (d) The results of the SEM. Arrows a and b indicate the Rayleigh waves and the converted waves caused by the free surface, respectively. The arrows c indicate the converted S-wave reflected from the inclined interface. Figure 6. Open in new tabDownload slide Snapshots of wave propagation in a two-layered model with a tilted interface at t = 1.0 s. (a), (b) The horizontal and vertical components of the displacement calculated with MFEM. (c), (d) The results of the SEM. Arrows a and b indicate the Rayleigh waves and the converted waves caused by the free surface, respectively. The arrows c indicate the converted S-wave reflected from the inclined interface. The seismograms recorded at the receivers are shown in figure 7. The fluctuation of the Rayleigh wave denoted by arrow a is observed at receiver R1. Arrows b and c indicate the head wave and the converted S-wave reflected from the inclined interface, respectively. The results of the MFEM are in good agreement with those of the SEM. The residuals are less than 1% and are thus negligible. However, the MFEM has an advantage over the SEM in computational efficiency (table 3). It is observed that although the CPU time required by the MFEM approximates to that by the SEM, the memory cost in the MFEM is only about 20% of that in the SEM. Figure 7. Open in new tabDownload slide Seismograms recorded at receivers R1 and R2 in the two-layered model. (a), (b) Horizontal and vertical displacement components at R1, respectively. (c), (d) Horizontal and vertical displacement components at R2, respectively. The blue and green lines represent the results of MFEM and SEM, respectively. The red lines indicate the residuals. Figure 7. Open in new tabDownload slide Seismograms recorded at receivers R1 and R2 in the two-layered model. (a), (b) Horizontal and vertical displacement components at R1, respectively. (c), (d) Horizontal and vertical displacement components at R2, respectively. The blue and green lines represent the results of MFEM and SEM, respectively. The red lines indicate the residuals. Table 3. Computational efficiency of MFEM and SEM in the second model. Methods . MFEM . SEM . Memory (KB) 341 100 1660 944 CPU time (s) 3875.97 3923.30 Methods . MFEM . SEM . Memory (KB) 341 100 1660 944 CPU time (s) 3875.97 3923.30 Open in new tab Table 3. Computational efficiency of MFEM and SEM in the second model. Methods . MFEM . SEM . Memory (KB) 341 100 1660 944 CPU time (s) 3875.97 3923.30 Methods . MFEM . SEM . Memory (KB) 341 100 1660 944 CPU time (s) 3875.97 3923.30 Open in new tab 3.3. A rugged surface model In the third numerical experiment, a model with a rugged free surface and interface is employed to validate the accuracy of the MFEM in complex media. Figure 8 shows the geometry of the model. The maximum height of the rugged free surface is about 350 m. The model is 2000 m wide and 1000 m deep. The thickness of the PML along the computational domain is 100 m. There are two irregular isotropic layers in the model. The upper layer is characterized by a P-wave velocity of 2000 m s-1, an S-wave velocity of 1250 m s-1 and a density of 2000 kg m-3. The lower layer is characterized by a P-wave velocity of 3000 m s-1, an S-wave velocity of 1578 m s-1 and a density of 2700 kg m-3. The source located at (1000 m, -30 m) is a Ricker wavelet with a central frequency of 15 Hz. Figure 8. Open in new tabDownload slide Geometry of the model with a rugged free surface and interface. The source indicated with a star is located at (1000 m, -30 m). Figure 8. Open in new tabDownload slide Geometry of the model with a rugged free surface and interface. The source indicated with a star is located at (1000 m, -30 m). Figure 9 shows snapshots of horizontal displacement at t = 0.3 s (figures 9(a), (c) and (e)) and t = 0.5 s (figures 9(b), (d) and (f)). The top, middle and bottom rows represent the wavefields calculated with the first-, second- and third-order MFEM, respectively. The wavefields become very complex due to the rugged free surface and interface. The arrows in (a) and (b) indicate the direct waves and the reflected waves from the free surface and interface. The arrows in (c) indicate the converted waves from the reflection and transmission of the interface. There is a clear numerical dispersion in the results of the first-order MFEM due to the poor accuracy. The dispersion is illustrated with arrows in figures 9(a) and (b). However, the numerical dispersion can be greatly reduced as the interpolation order increases. The resultant wavefields become approximately close to the physically reasonable results, which are shown in figures 9(c)–(f). Figure 9. Open in new tabDownload slide Snapshots of horizontal displacement simulated with the MFEM in different orders: (a), (c) and (e) are the results at t = 0.3 s and (b), (d) and (f) are the results at t = 0.5 s. The top, middle and bottom rows represent wavefields calculated with the first-, second- and third-order MFEM, respectively. Figure 9. Open in new tabDownload slide Snapshots of horizontal displacement simulated with the MFEM in different orders: (a), (c) and (e) are the results at t = 0.3 s and (b), (d) and (f) are the results at t = 0.5 s. The top, middle and bottom rows represent wavefields calculated with the first-, second- and third-order MFEM, respectively. The computational costs of the FEM and MFEM in the first-, second- and third-order are listed in table 4. In the simulation, all the programs are serial and run on the same personal computer. The MFEM needs less memory than the FEM. With the increase in interpolation order, the difference between the memory requirements of the FEM and MFEM becomes significant. The memory costs in the MFEM with first-, second- and third-order interpolation are 91%, 45% and 32% of those in the FEM, respectively. However, the MFEM solves the equations element-by-element and thus takes up more CPU time. The costs of CPU time in the first-, second- and third-order MFEM are about 2.09, 1.39 and 1.19 times as that in the FEM, respectively. Therefore, the computational efficiency of the MFEM becomes higher with the increase in numerical precision. Table 4. Computational efficiency of the first- to third-order FEM and MFEM in the third model. . FEM . MFEM . Methods . Memory (KB) . CPU time (s) . Memory (KB) . CPU time (s) . First-order 19 564 35.49 17 884 74.48 Second-order 91 672 233.12 41 576 326.24 Third-order 276 060 742.28 88 452 885.33 . FEM . MFEM . Methods . Memory (KB) . CPU time (s) . Memory (KB) . CPU time (s) . First-order 19 564 35.49 17 884 74.48 Second-order 91 672 233.12 41 576 326.24 Third-order 276 060 742.28 88 452 885.33 Open in new tab Table 4. Computational efficiency of the first- to third-order FEM and MFEM in the third model. . FEM . MFEM . Methods . Memory (KB) . CPU time (s) . Memory (KB) . CPU time (s) . First-order 19 564 35.49 17 884 74.48 Second-order 91 672 233.12 41 576 326.24 Third-order 276 060 742.28 88 452 885.33 . FEM . MFEM . Methods . Memory (KB) . CPU time (s) . Memory (KB) . CPU time (s) . First-order 19 564 35.49 17 884 74.48 Second-order 91 672 233.12 41 576 326.24 Third-order 276 060 742.28 88 452 885.33 Open in new tab 4. Conclusions In this study, a MFEM is developed for elastic wave modeling and has the following advantages. Firstly, the triangular elements are utilized and thus can fit the rugged free surface and discontinuous interfaces. Secondly, the kernel matrices of the global stiffness matrix, not the whole matrix, are stored, thus significantly reducing memory requirements. In addition, an element-by-element scheme is adopted, which avoids the assembly of large sparse matrices and makes parallel programming easier. Compared with the constant bandwidth storage scheme, the CSR methods of the MFEM save a large amount of computer memory, and the efficiency increases with the spatial interpolation order. The MFEM is feasible for simulating surface waves and all types of body waves accurately. Although a little more time is taken to calculate the matrix–vector product element-by-element, the cost is still acceptable. The results of wave propagation in a model with a rugged free surface suggest that the higher the order is, the better the performance of MFEM in saving both memory requirements and CPU times. MFEM can be a good choice for wave equation-based seismic imaging methods, owing to the higher computational efficiency, lower memory requirement and the flexibility in subdividing complex geometries. Acknowledgments The authors thank the four anonymous reviewers for their detailed suggestions, which greatly improved the quality of the paper. The C programs of 2D FEM and SEM are from Dr Shaolin Liu. This research was supported by the Natural Science Foundation of China (grant no. 41130418) and the National High Technology Research and Development Program (863 Program) of China (grant no. 2013AA064202). Appendix The empirical expressions for the damping coefficients, dx(x) and dz(z), are given by Collino and Tsogka (2001) dx(x)=-(n+1)⋅log(R)⋅VP2δ⋅xδ2dz(z)=-(n+1)⋅log(R)⋅VP2δ⋅zδ2 where δ is the thickness of the PML absorbing layer, x or z is the distance to the computational boundary in the PML region, n is a constant that controls the attenuation rate (which generally takes a value of 2 or 3) and R is the theoretical reflection coefficient which ranges from 10-6 to 10-2 in general. References Alford R M , Kelly K R , Boore D M . , 1974 Accuracy of finite-difference modeling of the acoustic wave equation , Geophysics , vol. 39 (pg. 834 - 842 ) 10.1190/1.1440470 Google Scholar Crossref Search ADS WorldCat Crossref Archer G C , Whalen T M . , 2005 Development of rotationally consistent diagonal mass matrices for plate and beam elements , Comput. Methods Appl. Mech. Eng. , vol. 194 (pg. 675 - 689 ) 10.1016/j.cma.2003.08.015 Google Scholar Crossref Search ADS WorldCat Crossref Bohlen T , Saenger E H . , 2006 Accuracy of heterogeneous staggered-grid finite-difference modeling of Rayleigh waves , Geophysics , vol. 71 (pg. T109 - T115 ) 10.1190/1.2213051 Google Scholar Crossref Search ADS WorldCat Crossref Chaljub E , Komatitsch D , Vilotte J P , Capdeville Y , Valette B , Festa G . , 2007 Spectral-element analysis in seismology , Adv. Geophys. , vol. 48 (pg. 365 - 419 ) 10.1016/S0065-2687(06)48007-9 Google Scholar Crossref Search ADS WorldCat Crossref Chen J B . , 2009 Lax–Wendroff and Nyström methods for seismic modeling , Geophys. Prospect. , vol. 57 (pg. 931 - 941 ) 10.1111/j.1365-2478.2009.00802.x Google Scholar Crossref Search ADS WorldCat Crossref Collino F , Tsogka C . , 2001 Application of the perfectly matched absorbing layer model to the linear elastodynamic problem in anisotropic heterogeneous media , Geophysics , vol. 66 (pg. 294 - 307 ) 10.1190/1.1444908 Google Scholar Crossref Search ADS WorldCat Crossref Cuthill E , McKee J . , 1969 Reducing the bandwidth of sparse symmetric matrices Proc. 1969 24th Natl Conf. of the ACM (pg. 157 - 172 ) De Hoop A T . , 1960 A modification of Cagniard’s method for solving seismic pulse problems , Appl. Sci. Res. , vol. 8 (pg. 349 - 356 ) 10.1007/BF02920068 Google Scholar Crossref Search ADS WorldCat Crossref Fornberg B . , 1988 The pseudospectral method: accurate representation of interfaces in elastic wave calculations , Geophysics , vol. 53 (pg. 625 - 637 ) 10.1190/1.1442497 Google Scholar Crossref Search ADS WorldCat Crossref Gibbs N E , Poole W G , Stockmeyer P K . , 1976 An algorithm for reducing the bandwidth and profile of a sparse matrix , SIAM J. Numer. Anal. , vol. 13 (pg. 236 - 250 ) Google Scholar Crossref Search ADS WorldCat Hesthaven J S , Warburton T . , 2008 , Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications New York Springer Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC He X , Yang D , Wu H . , 2015 A weighted Runge–Kutta discontinuous Galerkin method for wavefield modeling , Geophys. J. Int. , vol. 200 (pg. 1389 - 1410 ) 10.1093/gji/ggu487 Google Scholar Crossref Search ADS WorldCat Crossref Jiang L , Pang Z . , 1979 , The Method and the Theoretical Basis of the Finite Element Beijing People’s Education Press (in Chinese) Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Kelly K R , Ward R W , Treitel S , Alford R M . , 1976 Synthetic seismograms: a finite-difference approach , Geophysics , vol. 41 (pg. 2 - 27 ) 10.1190/1.1440605 Google Scholar Crossref Search ADS WorldCat Crossref Komatitsch D , Vilotte J P . , 1998 The spectral element method: an efficient tool to simulate the seismic response of 2D and 3D geological structures , Bull. Seismol. Soc. Am. , vol. 88 (pg. 368 - 392 ) OpenURL Placeholder Text WorldCat Komatitsch D , Tromp J . , 2002 Spectral-element simulations of global seismic wave propagation—I validation , Geophys. J. Int. , vol. 149 (pg. 390 - 412 ) 10.1046/j.1365-246X.2002.01653.x Google Scholar Crossref Search ADS WorldCat Crossref Komatitsch D , Tromp J . , 2003 A perfectly matched layer absorbing boundary condition for the second-order seismic wave equation , Geophys. J. Int. , vol. 154 (pg. 146 - 153 ) 10.1046/j.1365-246X.2003.01950.x Google Scholar Crossref Search ADS WorldCat Crossref Komatitsch D , Erlebacher G , Göddeke D , Michéa D . , 2010 High-order finite-element seismic wave propagation modeling with MPI on a large GPU cluster , J. Comput. Phys. , vol. 229 (pg. 7692 - 7714 ) 10.1016/j.jcp.2010.06.024 Google Scholar Crossref Search ADS WorldCat Crossref Li X F , Zhu T , Zhang M G , Long G H . , 2010 Seismic scalar wave equation with variable coefficients modeling by a new convolutional differentiator , Comput. Phys. Commun. , vol. 181 (pg. 1850 - 1858 ) 10.1016/j.cpc.2010.07.009 Google Scholar Crossref Search ADS WorldCat Crossref Li X F , Li Y Q , Zhang M G , Zhu T . , 2011 Scalar seismic-wave equation modeling by a multisymplectic discrete singular convolution differentiator method , Bull. Seismol. Soc. Am. , vol. 101 (pg. 1710 - 1718 ) 10.1785/0120100266 Google Scholar Crossref Search ADS WorldCat Crossref Li X F , Wang W S , Lu M W , Zhang M G , Li Y Q . , 2012 Structure-preserving modeling for elastic waves: a symplectic discrete singular convolution differentiator method , Geophys. J. Int. , vol. 188 (pg. 1382 - 1392 ) 10.1111/j.1365-246X.2011.05344.x Google Scholar Crossref Search ADS WorldCat Crossref Liu S L , Li X F , Wang W S , Liu Y S . , 2014a A mixed-grid finite element method with PML absorbing boundary conditions for seismic wave modeling , J. Geophys. Eng. , vol. 11 055009 10.1088/1742-2132/11/5/055009 OpenURL Placeholder Text WorldCat Crossref Liu S L , Li X F , Liu Y S , Zhu T , Zhang M G . , 2014b Dispersion analysis of triangle-based finite element method for acoustic and elastic wave simulation , Chin. J. Geophys. , vol. 57 (pg. 2620 - 2630 ) (in Chinese) OpenURL Placeholder Text WorldCat Liu S L , Li X F , Wang W S , Zhu T . , 2015a A symplectic RKN scheme for solving elastic wave equations , Chin. J. Geophys. , vol. 58 (pg. 1355 - 1366 ) (in Chinese) OpenURL Placeholder Text WorldCat Liu S L , Li X F , Wang W S , Zhu T . , 2015b Source wavefield reconstruction using a linear combination of the boundary wavefield in reverse time migration , Geophysics , vol. 80 (pg. S203 - S212 ) 10.1190/geo2015-0109.1 Google Scholar Crossref Search ADS WorldCat Crossref Liu Y S , Liu S L , Zhang M G , Ma D T . , 2012 An improved perfectly matched layer absorbing boundary condition for second order elastic wave equation , Prog. Geophys. , vol. 27 (pg. 2113 - 2122 ) OpenURL Placeholder Text WorldCat Liu Y S , Teng J W , Liu S L , Xu T . , 2013 Explicit finite element method with triangle meshes stored by sparse format and its perfectly matched layers absorbing boundary condition , Chin. J. Geophys. , vol. 56 (pg. 3085 - 3099 ) (in Chinese) OpenURL Placeholder Text WorldCat Liu Y , Teng J , Lan H , Si X , Ma X . , 2014 A comparative study of finite element and spectral element methods in seismic wavefield modeling , Geophysics , vol. 79 (pg. T91 - 104 ) 10.1190/geo2013-0018.1 Google Scholar Crossref Search ADS WorldCat Crossref Ma X , Yang D , Liu F . , 2011 A nearly analytic symplectically partitioned Runge–Kutta method for 2D seismic wave equations , Geophys. J. Int. , vol. 187 (pg. 480 - 496 ) 10.1111/j.1365-246X.2011.05158.x Google Scholar Crossref Search ADS WorldCat Crossref Marfurt K J . , 1984 Accuracy of finite-difference and finite-element modeling of the scalar and elastic wave equations , Geophysics , vol. 49 (pg. 533 - 549 ) 10.1190/1.1441689 Google Scholar Crossref Search ADS WorldCat Crossref Moczo P , Kristek J , Halada L . , 2000 3D fourth-order staggered-grid finite-difference schemes: stability and grid dispersion , Bull. Seismol. Soc. Am. , vol. 90 (pg. 587 - 603 ) 10.1785/0119990119 Google Scholar Crossref Search ADS WorldCat Crossref Padovani E , Priolo E , Seriani G . , 1994 Low and high order finite element method: experience in seismic modeling , J. Comput. Acoust. , vol. 2 (pg. 371 - 422 ) 10.1142/S0218396X94000233 Google Scholar Crossref Search ADS WorldCat Crossref Richter G R . , 1994 An explicit finite element method for the wave equation , Appl. Numer. Math. , vol. 16 (pg. 65 - 80 ) 10.1016/0168-9274(94)00048-4 Google Scholar Crossref Search ADS WorldCat Crossref Shi R , Wang S , Zhao J . , 2012 An unsplit complex-frequency-shifted PML based on matched Z-transform for FDTD modeling of seismic wave equations , J. Geophys. Eng. , vol. 9 (pg. 218 - 229 ) 10.1088/1742-2132/9/2/218 Google Scholar Crossref Search ADS WorldCat Crossref Tessmer E , Kosloff D . , 1994 3D elastic modeling with surface topography by a Chebychev spectral method , Geophysics , vol. 59 (pg. 464 - 473 ) 10.1190/1.1443608 Google Scholar Crossref Search ADS WorldCat Crossref Virieux J . , 1984 SH-wave propagation in heterogeneous media: velocity–stress finite-difference method , Geophysics , vol. 49 (pg. 1933 - 1942 ) 10.1190/1.1441605 Google Scholar Crossref Search ADS WorldCat Crossref Virieux J . , 1986 P-SV wave propagation in heterogeneous media: velocity-stress finite-difference method , Geophysics , vol. 51 (pg. 889 - 901 ) 10.1190/1.1442147 Google Scholar Crossref Search ADS WorldCat Crossref Wang Y , Liang W , Nashed Z , Li X , Liang G , Yang C . , 2014 Seismic modeling by optimizing regularized staggered-grid finite-difference operators using a time-space-domain dispersion-relationship-preserving method , Geophysics , vol. 79 (pg. T277 - T285 ) 10.1190/geo2014-0078.1 Google Scholar Crossref Search ADS WorldCat Crossref Wu S R . , 2006 Lumped mass matrix in explicit finite element method for transient dynamics of elasticity , Comput. Methods Appl. Mech. Eng. , vol. 195 (pg. 5983 - 5994 ) 10.1016/j.cma.2005.10.008 Google Scholar Crossref Search ADS WorldCat Crossref Yang D , Wang N , Chen S , Song G . , 2009 An explicit method based on the implicit Runge–Kutta algorithm for solving wave equations , Bull. Seismol. Soc. Am. , vol. 99 (pg. 3340 - 3354 ) 10.1785/0120080346 Google Scholar Crossref Search ADS WorldCat Crossref Zhang M G , Wang M Y , Li X F , Yang X C , Wang L . , 2002 Finite element forward modeling of anisotropic elastic waves , Prog. Geophys. , vol. 17 (pg. 384 - 389 ) (in Chinese) OpenURL Placeholder Text WorldCat Zhang W , Shen Y . , 2010 Unsplit complex frequency shifted PML implementation using auxiliary differential equation for seismic wave modeling , Geophysics , vol. 75 (pg. T141 - T154 ) 10.1190/1.3463431 Google Scholar Crossref Search ADS WorldCat Crossref Zhang W , Zhang Z , Chen X . , 2012 Three-dimensional elastic wave numerical modeling in the presence of surface topography by a collocated-grid finite-difference method on curvilinear grids , Geophys. J. Int. , vol. 190 (pg. 358 - 378 ) 10.1111/j.1365-246X.2012.05472.x Google Scholar Crossref Search ADS WorldCat Crossref © 2017 Sinopec Geophysical Research Institute
TI - Seismic wavefield simulation by a modified finite element method with a perfectly matched layer absorbing boundary
JF - Journal of Geophysics and Engineering
DO - 10.1088/1742-2140/aa6b31
DA - 2017-08-01
UR - https://www.deepdyve.com/lp/oxford-university-press/seismic-wavefield-simulation-by-a-modified-finite-element-method-with-o85b1Y3K36
SP - 852
VL - 14
IS - 4
DP - DeepDyve
ER -