The degree of a complete polynomial
Abstract
In this chapter, a new efficient highorder finite volume method for 3D elastic modelling on unstructured meshes is developed. The stencil for the highorder polynomial reconstruction is generated by subdividing the relative coarse tetrahedrons. The reconstruction on the stencil is performed by using cellaveraged quantities represented by the hierarchical orthonormal basis functions. Unlike the traditional highorder finite volume method, the new method has a very local property like the discontinuous Galerkin method. Furthermore, it can be written as an innersplit computational scheme which is beneficial to reducing computational amount. The reconstruction matrix is invertible and remains unchanged for all tetrahedrons, and thus it can be precomputed and stored before time evolution. These special advantages facilitate the parallelization and highorder computations. The highorder accuracy in time is obtained by the RungeKutta method. Numerical computations including a 3D real model with complex topography demonstrate the effectiveness and good adaptability to complex topography.
Keywords
 numerical solutions
 computational seismology
 3D elastic wave
 wave propagation
 highorder finite volume method
 unstructured meshes
1. Introduction
Wave propagation based on wave equations has important applications in geophysics. It is usually used as a powerful tool to detect the structures of reservoir. Thus solving wave equations efficiently and accurately is always an important research topic. There are several types of numerical methods to solve wave equations, for example, the finite difference (FD) method [1, 2], the pseudospectral (PS) method [3, 4], the finite element (FE) method [5, 6, 7, 8, 9], the spectral element (SE) method [10, 11, 12, 13, 14], the discontinuous Galerkin (DG) method [15, 16, 17, 18], and the finite volume (FV) method [19, 20, 21, 22]. Each numerical method has its own inherent advantages and disadvantages. For example, the FD method is efficient and relatively easy to implement, but the inherent restriction of using regular meshes limits its application to complex topography. The FE method has good adaptability to complex topography, but it has huge computational cost. In this chapter, the FV method is the key consideration.
In order to simulate wave propagation on unstructured meshes efficiently, the FV method is a good choice due to its high computational efficiency and good adaptability to complex geometry. In this chapter an efficient FV method for 3D elastic wave simulation on unstructured meshes is developed. It incorporates some nice features from the DG and FV methods [15, 16, 17, 19, 20, 23] and the spectral FV (SFV) method [24, 25, 26]. In our method, the computational domain is first meshed with relative coarse tetrahedral elements in 3D or triangle elements in 2D. Then, each element is further divided as a collection of finer subelements to form a stencil. The highorder polynomial reconstruction is performed on this stencil by using local cellaveraged values on the finer elements. The resulting reconstruction matrix on all coarse elements remains unchanged, and it can be precomputed before time evolution. Moreover, the method can be written as an innersplit computational scheme. These two advantages of our method are very beneficial to enhancing the parallelization and reducing computational cost.
The rest of this chapter is organized as follows. In Section 2, the theory is described in detail. In Section 3, numerical results are given to illustrate the effectiveness of our method. Finally, the conclusion is given in Section 4.
2. Theory
2.1 The governing equation
The threedimensional (3D) elastic wave equation with external sources in velocitystress formulation can be written as the following system [1, 15]:
where
where
The propagation velocities of the elastic waves are determined by the eigenvalues
where
are the velocities of the compression (
2.2 The generation of a stencil
Suppose that the 3D computational domain
In practical computations, the integrals in the FV scheme on physical tetrahedral element
and its corresponding inverse transformation by
The detailed expressions of the transformations (6) and (7) will be given in Section 2.5.
Inside each
where
In order to construct a highorder polynomial, we need to choose a stencil. Traditionally, the elements being adjacent to the element
Let
In our algorithm, we guarantee

1  2  3  4 


4  10  20  35 

2  3  3  4 

8  27  27  64 
2.3 The highorder polynomial reconstruction
The highorder polynomial is reconstructed in each element
where
to reconstruct a highorder polynomial, where
To solve the reconstruction problem, inspired by the DG method [15, 16, 17, 23, 28, 29], we use hierarchical orthogonal basis functions. The basis functions
Transforming equation (12) in the physical coordinate system
where
The integration in Eq. (14) over
where
and
We need at least
From the orthogonality of basis functions and the property of Eq. (13), we remark that Eq. (15) is subject to the following constraint condition [27]:
With the constraint, Eq. (15) is solved by the Lagrange multiplier method [19, 20, 27]. And the system can be written as
where
The coefficient matrix on the lefthand side of Eq. (19) is the socalled reconstruction matrix [19, 20].
2.4 The spatial discrete formulation
We now derive the semidiscrete finite volume scheme based on Eqs. (2) and (8). Integrating over each subelement
Using Eq. (8) and integration by parts yield
where
where
where
And
where
Inserting Eqs. (23) into (22) and rewriting the result into a splitting form of easy computation in the reference system
with
and
where
where

1  2  3 









2.5 The time discretization
Equation (27) is in fact a semidiscrete ordinary differential equation (ODE) system. In order to solve it formally, we denote the spatial semidiscrete part in Eq. (27) by a linear operator
Traditionally, the classic fourthorder explicit RK (ERK) method
can be applied to advance
As we can see the LSERK only requires one additional storage level, while ERK has four. The coefficients required in Eq. (34) are listed in Table 3 [30].





1  0  0.1496590219992291  0 
2  −0.4178904744998519  0.3792103129996273  0.1496590219992291 
3  −1.1921516946426769  0.8229550293869817  0.3704009573642048 
4  −1.6977846924715279  0.6994504559491221  0.6222557631344432 
5  −1.5141834442571558  0.1530572479681520  0.9582821306746903 
As to the stability condition, it is controlled by the CourantFriedrichsLewy (CFL) condition [15, 19];
where
The absorbing boundary conditions (ABCs) in computations are required as the computational domain is finite. There are two typical ABCs to be adopted here. One is flux type ABCs [16, 19]. That is to say, the following numerical flux in Eq. (23) at all tetrahedral faces that coincide with domain boundary
which allows only for outgoing waves and is equivalent to the first order ABCs. Though the absorbing effects of this method vary the angles of incidence, it is still effective in many cases [19]. The advantage of this type ABCs is that it merged into the FVM framework naturally and there is almost no additional computational cost. Another type is the perfectly matched layer (PML) technique originally developed by [31], which is very popular in recent more 10 years.
2.6 Coordinate transformation
The transformation between different coordinate systems is frequently used. For ease of reading, we present the formulations here. Let
then the transformation from
where
Note that
The coordinate transformation from the second reference coordinate
then the transform from
which is the determinant of the Jacobian matrix of the transformation being equal to six times the volume of the subelement
3. Numerical computations
In this section we give three numerical examples to illustrate the performance of the developed method above. The convergence test of the proposed method can be found in [27]. Though the method is developed for the 3D case, it can be simplified to 2D without essential difficulty. The principle is the same. The first example is a test for a 2D model with uneven topography. The other two examples are for two 3D models.
Example 1. The first example is a twolayered model with the inclined interface shown in Figure 3a. The range of the model is
where
is applied, where
Example 2. The second example is a cuboid model. The physical size of the model is
It is applied to the
Example 3. The third example is a real geological model in China. As shown in Figure 8a, it has a very complex topography. The physical scope of the model is
4. Conclusions
A new efficient highorder finite volume method for the 3D elastic wave simulation on unstructured meshes has been developed. It combines the advantages of the DG method and the traditional FV method. It adapts irregular topography very well. The reconstruction stencil is generated by refining each coarse tetrahedron which can be implemented effectively for all tetrahedrons whether they are internal or boundary elements. The hierarchical orthogonal basis functions are exploited to perform the highorder polynomial reconstruction on the stencil. The resulting reconstruction matrix remains unchanged for all tetrahedrons and can be precomputed and stored before time evolution. The method preserves a very local property like the DG method, while it has high computational efficiency like the FV method. These advantages facilitate 3D largescale parallel computations. Numerical computations including a 3D real physical model show its good performance. The method also can be expected to solve other linear hyperbolic equations without essential difficulty.
Acknowledgments
I appreciate Dr. Y. Zhuang, Prof. Chung, and Dr. L. Zhang very much for their important help and cooperation. This work is supported by the National Natural Science Foundation of China under the grant number 11471328 and 51739007. It is also partially supported by the National Center for Mathematics and Interdisciplinary Sciences, Chinese Academy of Sciences.
References
 1.
Minkoff SE. Spatial parallelism of a 3D finite difference velocitystress elastic wave propagation code. SIAM Journal on Scientific Computing. 2002; 24 :119. DOI: 10.1137/S1064827501390960  2.
Virieux J. PSV wave propagation in heterogeneous media: Velocitystress finitedifference method. Geophysics. 1986; 51 :889901. DOI: 10.1190/1.1442147  3.
Klin P, Priolo E, Seriani G. Numerical simulation of seismic wave propagation in realistic 3D geomodels with a Fourier pseudo spectral method. Geophysical Journal International. 2010; 183 :905922. DOI: 10.1111/j.1365246X.2010.04763.x  4.
Zhang W. Stability conditions for wave simulation in 3D anisotropic media with the pseudospectral method. Communications in Computational Physics. 2012; 12 :703720. DOI: 10.4208/cicp.120610.090911a  5.
Bécache E, Joly P, Tsogka C. A new family of mixed finite elements for the linear elastodynamic problem. SIAM Journal on Numerical Analysis. 2002; 39 :21092132. DOI: 10.2307/4101053  6.
Cohen G, Joly P, Roberts JE, Tordjman N. Higher order triangular finite elements with mass lumping for the wave equations. SIAM Journal on Numerical Analysis. 2001; 38 :20472078. DOI: 10.1137/s0036142997329554  7.
Cohen G, Fauqueux S. Mixed finite elements with masslumping for the transient wave equation. Journal of Computational Acoustics. 2000; 8 :171188. DOI: 10.1142/s0218396x0000011x  8.
Cohen G, Fauqueux S. Mixed spectral finite elements for the linear elasticity system in unbounded domains. SIAM Journal on Scientific Computing. 2005; 26 :864884. DOI: 10.1137/S1064827502407457  9.
Zhang W, Chung E, Wang C. Stability for imposing absorbing boundary conditions in the finite element simulation of acoustic wave propagation. Journal of Computational Mathematics. 2014; 32 :120. DOI: 10.4208/jcm.1310m3942  10.
Dubiner M. Spectral methods on triangles and other domains. Journal of Scientific Computing. 1991; 6 :345390. DOI: 10.1007/BF01060030  11.
Komatitsch D, Tromp J. Introduction to the spectral element method for threedimensional seismic wave propagation. Geophysical Journal International. 1999; 139 :806822. DOI: 10.1046/j.1365246x.1999.00967.x  12.
Komatitsch D, Martin R, Tromp J, Taylor MA, Wingate BA. Wave propagation in 2D elastic media using a spectral element method with triangles and quadrangles. Journal of Computational Acoustics. 2001; 9 :703718. DOI: 10.1142/S0218396X01000796  13.
Komatitsch D, Tromp J. Spectralelement simulations of global seismic wave propagation—I. Validation. Geophysical Journal International. 2002; 149 :390412. DOI: 10.1046/j.1365246X.2002.01653.x  14.
Seriani G. 3D largescale wave propagation modeling by spectralelement method on Cray T3E multiprocessor. Computer Methods in Applied Mechanics and Engineering. 1998; 164 :235247. DOI: 10.1016/S00457825(98)000577  15.
Dumbser M, Käser M. An arbitrary highorder discontinuous Galerkin method for elastic waves on unstructured meshesII: The threedimensional isotropic case. Geophysical Journal International. 2006; 167 :319336. DOI: 10.1111/j.1365246X.2006.03120.x  16.
Käser M, Dumbser M. An arbitrary highorder discontinuous Galerkin method for elastic waves on unstructured meshes. I: The twodimensional isotropic case with external source terms. Geophysical Journal International. 2006; 166 :855877. DOI: 10.1111/j.1365246X.2006.03051.x  17.
Käser M, Dumbser M, Puente J, Igel H. An arbitrary highorder discontinuous Galerkin method for elastic waves on unstructured meshes—III: Viscoelastic attenuation. Geophysical Journal International. 2007; 168 :224242. DOI: 10.1111/j.1365246X.2006.03193.x  18.
Ye R, de Hoop MV, Petrovitch CL, PyrakNolte LJ, Wilcox LC. A discontinuous Galerkin method with a modified penalty flux for the propagation and scattering of acoustoelastic waves. Geophysical Journal International. 2016; 205 :12671289. DOI: 10.1093/gji/ggw070  19.
Dumbser M, Käser M, de la Puente J. Arbitrary highorder finite volume schemes for seismic wave propagation on unstructured meshes in 2D and 3D. Geophysical Journal International. 2007; 171 :665694. DOI: 10.1111/j.1365246X.2007.03421.x  20.
Dumbser M, Käser M. Arbitrary high order nonoscillatory finite volume schemes on unstructured meshes for linear hyperbolic systems. Journal of Computational Physics. 2007; 221 :693723. DOI: 10.1016/j.jcp.2006.06.043  21.
Käser M, Iske A. ADER schemes on adaptive triangular meshes for scalar conservation laws. Journal of Computational Physics. 2005; 205 :486508. DOI: 10.1016/j.jcp.2004.11.015  22.
Zhang W, Zhuang Y, Chung ET. A new spectral finite volume method for elastic modelling on unstructured meshes. Geophysical Journal International. 2016; 206 :292307. DOI: 10.1093/gji/ggw148  23.
Dumbser M, Käser M, Toro EF. An arbitrary highorder discontinuous Galerkin method for elastic waves on unstructured meshes—V: Local time stepping and padaptivity. Geophysical Journal International. 2007; 171 :695717. DOI: 10.1111/j.1365246X.2007.03427.x  24.
Liu Y, Vinokur M, Wang ZJ. Spectral (finite) volume method for conservation laws on unstructured grids. V: Extension to threedimensional systems. Journal of Computational Physics. 2006; 212 :454472. DOI: 10.1016/j.jcp.2003.09.012  25.
Wang ZJ. Spectral (finite) volume method for conservation laws on unstructured grids: Basic formulation. Journal of Computational Physics. 2002; 178 :210251. DOI: 10.1006/jcph.2002.7041  26.
Wang ZJ, Liu Y. Spectral (finite) volume method for conservation laws on unstructured grids. II: Extension to twodimensional scalar equation. Journal of Computational Physics. 2002; 179 :665697. DOI: 10.1006/jcph.2002.7082  27.
Zhang W, Zhuang Y, Zhang L. A new highorder finite volume method for 3D elastic wave simulation on unstructured meshes. Journal of Computational Physics. 2017; 340 :534555. DOI: 10.1016/j.jcp.2017.03.050  28.
Puente J, Käser M, Dumbser M, Igel H. An arbitrary highorder discontinuous Galerkin method for elastic waves on unstructured meshesIV: Anisotropy. Geophysical Journal International. 2007; 169 :12101228. DOI: 10.1111/j.1365246X.2007.03381.x  29.
Puente J, Dumbser M, Käser M, Igel H. Discontinuous Galerkin method for propagation in poroelastic media. Geophysics. 2008; 73 :T77T97. DOI: 10.1190/1.2965027  30.
Hesthaven JS, Warburton T. Nodal Discontinuous Galerkin Methods. New York: SpringerVerlag; 2008. 502 p. DOI: 10.1007/9780387720678  31.
Bérenger JP. A perfectly matched layer for the absorption of electromagnetic waves. Journal of Computational Physics. 1994; 114 :185200. DOI: 10.1006/jcph.1996.0181  32.
Collino F, Tsogka C. Application of the perfectly matched absorbing layer model to the linear elastodynamic problem in anisotropic heterogeneous media. Geophysics. 2001; 66 :294307. DOI: 10.1190/1.1444908