A High-Order Finite Volume Method for 3D Elastic Modelling on Unstructured Meshes

Abstract

In this chapter, a new efficient high-order finite volume method for 3D elastic modelling on unstructured meshes is developed. The stencil for the high-order polynomial reconstruction is generated by subdividing the relative coarse tetrahedrons. The reconstruction on the stencil is performed by using cell-averaged quantities represented by the hierarchical orthonormal basis functions. Unlike the traditional high-order finite volume method, the new method has a very local property like the discontinuous Galerkin method. Furthermore, it can be written as an inner-split 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 pre-computed and stored before time evolution. These special advantages facilitate the parallelization and high-order computations. The high-order accuracy in time is obtained by the Runge-Kutta method. Numerical computations including a 3D real model with complex topography demonstrate the effectiveness and good adaptability to complex topography.

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 pseudo-spectral (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 high-order polynomial reconstruction is performed on this stencil by using local cell-averaged values on the finer elements. The resulting reconstruction matrix on all coarse elements remains unchanged, and it can be pre-computed before time evolution. Moreover, the method can be written as an inner-split 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 three-dimensional (3D) elastic wave equation with external sources in velocity-stress formulation can be written as the following system [1, 15]:

$$\left\{\begin{array}{l}\frac{\partial {\sigma }_{\mathit{xx}}}{\partial t}-\left(\lambda +2\mu \right)\frac{\partial u}{\partial x}-\lambda \frac{\partial v}{\partial y}-\lambda \frac{\partial w}{\partial z}=g_1,\\ \frac{\partial {\sigma }_{\mathit{yy}}}{\partial t}-\lambda \frac{\partial u}{\partial x}-\left(\lambda +2\mu \right)\frac{\partial v}{\partial y}-\lambda \frac{\partial w}{\partial z}=g_2,\\ \frac{\partial {\sigma }_{\mathit{zz}}}{\partial t}-\lambda \frac{\partial u}{\partial x}-\lambda \frac{\partial v}{\partial y}-\left(\lambda +2\mu \right)\frac{\partial w}{\partial z}=g_3,\\ \frac{\partial {\sigma }_{\mathit{xy}}}{\partial t}-\mu \left(\frac{\partial v}{\partial x}+\frac{\partial u}{\partial y}\right)=g_4,\\ \frac{\partial {\sigma }_{\mathit{yz}}}{\partial t}-\mu \left(\frac{\partial v}{\partial z}+\frac{\partial w}{\partial y}\right)=g_5,\\ \frac{\partial {\sigma }_{\mathit{xz}}}{\partial t}-\mu \left(\frac{\partial u}{\partial z}+\frac{\partial w}{\partial x}\right)=g_6,\\ \rho \frac{\partial u}{\partial t}-\frac{\partial {\sigma }_{\mathit{xx}}}{\partial x}-\frac{\partial {\sigma }_{\mathit{xy}}}{\partial y}-\frac{\partial {\sigma }_{\mathit{xz}}}{\partial z}=\rho g_7,\\ \rho \frac{\partial v}{\partial t}-\frac{\partial {\sigma }_{\mathit{xy}}}{\partial x}-\frac{\partial {\sigma }_{\mathit{yy}}}{\partial y}-\frac{\partial {\sigma }_{\mathit{yz}}}{\partial z}=\rho g_8,\\ \rho \frac{\partial w}{\partial t}-\frac{\partial {\sigma }_{\mathit{xz}}}{\partial x}-\frac{\partial {\sigma }_{\mathit{yz}}}{\partial y}-\frac{\partial {\sigma }_{\mathit{zz}}}{\partial z}=\rho g_9,\end{array}\right.$$ E1

where $u$ , $v$ , and $w$ are the wavefield of particle velocities in $x$ , $y$ , and $z$ directions, respectively; $\lambda$ and $\mu$ are the Lamé coefficients and $\rho$ is the density; $g_i\left(x,y,z,t\right)$ are the known sources; ${\sigma }_{\mathit{xx}}$ , ${\sigma }_{\mathit{yy}}$ , and ${\sigma }_{\mathit{zz}}$ are the normal stress components while ${\sigma }_{\mathit{xy}}$ , ${\sigma }_{\mathit{xz}}$ , and ${\sigma }_{\mathit{yz}}$ are the shear stresses. For the convenient of discussion, we rewrite Eq. (1) as the following compact form:

$$\frac{\partial \mathbf{u}}{\partial t}+A\frac{\partial \mathbf{u}}{\partial x}+B\frac{\partial \mathbf{u}}{\partial y}+C\frac{\partial \mathbf{u}}{\partial z}=\mathbf{g},$$ E2

where $\mathbf{g}={\left(g_1,\cdots ,g_9\right)}^T$ , $\mathbf{u}={\left({\sigma }_{\mathit{xx}},{\sigma }_{\mathit{yy}},{\sigma }_{\mathit{zz}},{\sigma }_{\mathit{xy}},{\sigma }_{\mathit{yz}},{\sigma }_{\mathit{xz}},u,v,w\right)}^T$ , and the matrices $A$ , $B$ , and $C$ are all $9\times 9$ matrices and can be obtained obviously [27].

The propagation velocities of the elastic waves are determined by the eigenvalues $s_i$ of matrices $A$ , $B$ , and $C$ and are given by

$$s_1=-v_p,s_2=-v_s,s_3=-v_s,s_4=s_5=s_6=0,s_7=v_s,s_8=v_s,s_9=v_p,$$ E3

where

$$v_p=\sqrt{\frac{\lambda +2\mu }{\rho }},v_s=\sqrt{\frac{\mu }{\rho }}$$ E4

are the velocities of the compression ( $P$ ) wave and the shear ( $S$ ) wave velocities, respectively.

2.2 The generation of a stencil

Suppose that the 3D computational domain $\mathrm{\Omega }$ is meshed by $N_E$ conforming tetrahedral elements $T^{\left(m\right)}$ :

$$\mathrm{\Omega }=\bigcup _{m=1}^{N_E}{T}^{\left(m\right)}.$$ E5

In practical computations, the integrals in the FV scheme on physical tetrahedral element $T^{\left(m\right)}$ are usually changed to be computed on its reference element. Figure 1 shows a physical tetrahedron $T^{\left(m\right)}$ in the physical system, and $x-y-z$ is transformed into a reference element $T_E$ in the reference system $\xi -\eta -\zeta$ . Let $\left(x_i,y_i,z_i\right)$ for $i=\mathrm{1,2,3,4}$ be the coordinates of physical element $T^{\left(m\right)}$ . The transformations between $x-y-z$ system and $\xi -\eta -\zeta$ system will be given in the final subsection of Section 2. For convenience, let $\mathbf{x}=\left(x,y,z\right)$ and $\mathbf{\xi }=\left(\xi ,\eta ,\zeta \right)$ . And denote the transformation from $\xi -\eta -\zeta$ system to $x-y-z$ system by

Figure 1.

The physical element $T^{\left(m\right)}$ (left) in the physical coordinate system $x-y-z$ is transformed into a reference element $T_E$ (right) in the reference coordinate system $\xi -\eta -\zeta$ .

$$\mathbf{x}=\mathbf{x}\left(T^{\left(m\right)},\mathbf{\xi }\right),$$ E6

and its corresponding inverse transformation by

$$\mathbf{\xi }=\mathbf{\xi }\left(T^{\left(m\right)},\mathbf{x}\right).$$ E7

The detailed expressions of the transformations (6) and (7) will be given in Section 2.5.

Inside each $T_E$ the solutions of Eq. (2) are approximated numerically by using a linear combination of polynomial basis functions ${\varphi }_l\left(\xi ,\eta ,\zeta \right)$ and the time-dependent coefficients ${\widehat{\mathbf{w}}}_l^{\left(m\right)}\left(t\right)$ :

$${\mathbf{u}}^{\left(m\right)}\left(\xi ,\eta ,\zeta ,t\right)=\sum _{l=1}^{N_p}{\widehat{\mathbf{w}}}_l^{\left(m\right)}\left(t\right){\varphi }_l\left(\xi ,\eta ,\zeta \right),$$ E8

where $N_p$ is the degree of freedom of a complete polynomial.

In order to construct a high-order polynomial, we need to choose a stencil. Traditionally, the elements being adjacent to the element $T^{\left(m\right)}$ are selected to form a stencil. In [20] three types of stencils, i.e., the central stencil, the primary sector stencil, and the reverse stencil, are investigated. These stencils usually choose $2N$ neighbors for the 3D reconstruction. Here $N$ is the degree of a complete polynomial. Due to geometrical issues, the reconstruction matrix resulting from these stencils may be not invertible. This may happen when all elements are aligned in a straight line [20]. In the following, we propose to partition $T^{\left(m\right)}$ or in fact its corresponding reference element $T_E$ into finer subelements to form a stencil. The subdivision algorithm guarantees the number of subelements is greater than the degrees of freedom of a complete polynomial. Moreover, this algorithm is easy to implement especially in 3D and for all elements whether they are internal or boundary elements.

Let $N_e$ be the number of subelements in $T^{\left(m\right)}$ after subdividing. For a complete polynomial of degree $N$ in 3D, a reconstruction requires at least $N_p$ subelements, where

$$N_p=\left(N+1\right)\left(N+2\right)\left(N+3\right)/6.$$ E9

In our algorithm, we guarantee $N_e$ is always greater than $N_p$ . As shown in Figure 2, we divide each edge of the reference element $T_E$ into $M$ uniform segments. Thus we have $N_e≔M^3$ tetrahedral subelements in $T_E$ . Note that a small subcubic in $T_E$ consists of six tetrahedrons. With the transformations of Eqs. (6) and (7), we denote all subelements in $T^{\left(m\right)}$ for a fixed $m$ by $T^{\left(m\left(k\right)\right)}$ for $k=1,\cdots ,N_e$ . In Table 1, the degree of a complete polynomial $N$ and its corresponding degrees of freedom $N_p$ are listed. Correspondingly, the number of $M$ and $N_e$ are also listed in Table 1. This algorithm for generating the stencil is easily implemented for all coarse tetrahedrons. Moreover, the reconstruction matrix resulting from this stencil is always invertible and remains unchanged for all elements $T^{\left(m\right)}$ for $m=1,\cdots ,N_E$ . Note that the reconstruction matrix may be not invertible if all elements are aligned on a straight line [15]. However, this will not happen here for our algorithm.

Figure 2.

The stencil obtained by subdividing the reference element $T_E$ into $M^{\mathit{3}}={\mathit{3}}^{\mathit{3}}$ tetrahedral subelements, where $M=\mathit{3}$ is the number of uniform segments on each edge of $T_E$ . Note that a small subcubic (red) in $T_E$ consists of six tetrahedrons.

$N$	1	2	3	4
$N_p$	4	10	20	35
$M$	2	3	3	4
$N_e$	8	27	27	64

Table 1.

The degree of a complete polynomial $N$ and its corresponding degrees of freedom $N_p$ are listed. Correspondingly, the number of uniform segments $M$ on each edge and the number of subelements $N_e$ are also listed.

2.3 The high-order polynomial reconstruction

The high-order polynomial is reconstructed in each element $T^{\left(m\right)}$ or $T_E$ . For the stencil designed above, we have

$$T^{\left(m\right)}=\bigcup _{k=1}^{N_e}{T}^{\left(m\left(k\right)\right)},$$ E10

where $k=1,\cdots ,N_e$ is the index for subelements in $T^{\left(m\right)}$ . The FV method will use the cell-averaged quantities, i.e.,

$${\overline{\mathbf{u}}}^{\left(m\left(k\right)\right)}=\frac{1}{\mid T^{\left(m\left(k\right)\right)}\mid }{\int }_{T^{\left(m\left(k\right)\right)}}{\mathbf{u}}^{\left(m\right)}\left(\mathbf{x}\right)\mathit{dV},k=1,\cdots ,N_e,$$ E11

to reconstruct a high-order polynomial, where $\mid T^{\left(m\left(k\right)\right)}\mid$ represents the volume of the subelement $T^{\left(m\left(k\right)\right)}$ . The time variable $t$ in ${\mathbf{u}}^{\left(m\right)}$ is omitted for discussion convenience. The reconstruction requires integral conservation for ${\mathbf{u}}^{\left(m\right)}$ in each subelement $T^{\left(m\left(k\right)\right)}$ , i.e.,

$$\begin{array}{r}{\int }_{T^{\left(m\left(k\right)\right)}}{\mathbf{u}}^{\left(m\right)}\left(\mathbf{x}\left(T^{\left(m\right)},\mathbf{\xi }\right)\right)\mathit{dV}=\mid T^{\left(m\left(k\right)\right)}\mid {\overline{\mathbf{u}}}^{\left(m\left(k\right)\right)},\\ \forall T^{\left(m\left(k\right)\right)}\subset T^{\left(m\right)},k=1,\cdots ,N_e.\end{array}$$ E12

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 ${\varphi }_l\left(\xi ,\eta ,\zeta \right)$ of a complete polynomial of degree $N$ ( $N=\mathrm{1,2,3,4}$ ) in the reference coordinate system can be found in [27]. We remark that the basis functions are orthonormal and satisfy the following property:

$${\int }_{T_E}{\varphi }_l\left(\xi ,\eta ,\zeta \right)\mathit{\text{dξdηdζ}}=\left\{\begin{array}{ll}\frac{\sqrt{6}}{6},& l=1,\\ 0,& l\ne 1.\end{array}\right.$$ E13

Transforming equation (12) in the physical coordinate system $x-y-z$ into the reference coordinate system $\xi -\eta -\zeta$ and noticing Eq. (8), we obtain

$$\begin{array}{r}\sum _{l=1}^{N_p}\left({\int }_{{\tilde{T}}^{\left(m\left(k\right)\right)}}{\varphi }_l\left(\xi ,\eta ,\zeta \right)\mathit{\text{dξdηdζ}}\right){\widehat{\mathbf{w}}}_l^{\left(m\right)}=\mid {\tilde{T}}^{\left(m\left(k\right)\right)}\mid {\overline{\mathbf{u}}}^{\left(m\left(k\right)\right)},\\ \forall {\tilde{T}}^{\left(m\left(k\right)\right)}\subset {\tilde{T}}^{\left(m\right)}=T_E,k=1,\cdots ,N_e,\end{array}$$ E14

where ${\tilde{T}}^{\left(m\right)}$ is in fact the reference element $T_E$ and ${\tilde{T}}^{\left(m\left(k\right)\right)}$ is the transformed element corresponding to the subelement $T^{\left(m\left(k\right)\right)}$ .

The integration in Eq. (14) over ${\tilde{T}}^{\left(m\left(k\right)\right)}$ in $\mathbf{\xi }$ system can be computed efficiently if it is performed over its reference element in a second reference system $\tilde{\mathbf{\xi }}$ . Denote the transformation from $\tilde{\mathbf{\xi }}$ to $\mathbf{\xi }$ and its inverse by $\mathbf{\xi }=\mathbf{\xi }\left({\tilde{T}}^{\left(m\left(k\right)\right)},\tilde{\mathbf{\xi }}\right)$ and $\tilde{\mathbf{\xi }}=\tilde{\mathbf{\xi }}\left({\tilde{T}}^{\left(m\left(k\right)\right)},\mathbf{\xi }\right)$ , respectively. Transforming Eq. (14) into $\tilde{\mathbf{\xi }}$ system and rewriting the result as a compact form, we have

$$G\widehat{\mathbf{w}}=\overline{\mathbf{u}},$$ E15

where $G$ is the $N_e\times N_p$ matrix with entries $G_{\mathit{kl}}$ given by

$$G_{\mathit{kl}}=\frac{1}{\mid T_E\mid }\left({\int }_{T_E}{\varphi }_l\left(\mathbf{\xi }\left({{\displaystyle \tilde{T}}}^{\left(m\left(k\right)\right)},{\displaystyle \tilde{\mathbf{\xi }}}\right)\right)d{\displaystyle \tilde{\xi }d\tilde{\eta }}d{\displaystyle \tilde{\zeta }}\right),k=1,\cdots ,N_e;l=1,\cdots ,N_p,$$ E16

and

$$\overline{\mathbf{u}}≔{\left({\overline{\mathbf{u}}}^{\left(m\left(1\right)\right)},{\overline{\mathbf{u}}}^{\left(m\left(2\right)\right)},\cdots ,{\overline{\mathbf{u}}}^{\left(m\left(N_e\right)\right)}\right)}^T,\widehat{\mathbf{w}}≔{\left({\widehat{\mathbf{w}}}_1^{\left(m\right)},{\widehat{\mathbf{w}}}_2^{\left(m\right)},\cdots ,{\widehat{\mathbf{w}}}_{N_p}^{\left(m\right)}\right)}^T.$$ E17

We need at least $N_p$ subelements in the stencil since the reconstructed number of degrees of freedom is $N_p$ . As listed in Table 1, $N_e$ subelements are used to form the stencil. Note that $N_e$ is definitely larger than $N_p$ , which is helpful to improve the reconstruction robustness [20, 21]. Thus Eq. (15) is an overdetermined problem. We use the constrained least squared technique to solve it.

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]:

$$\sqrt{6}{\widehat{\mathbf{w}}}_1^{\left(m\right)}=\sum _{k=1}^{N_e}\frac{{\overline{\mathbf{u}}}^{\left(m\left(k\right)\right)}}{N_e}.$$ E18

With the constraint, Eq. (15) is solved by the Lagrange multiplier method [19, 20, 27]. And the system can be written as

$$\left(\begin{array}{ll}2G^{T}G& -R^T\\ R& 0\end{array}\right)\left(\begin{array}{l}\widehat{\mathbf{w}}\\ {\lambda }_p\end{array}\right)=\left(\begin{array}{l}2G^T\overline{\mathbf{u}}\\ \tilde{R}\overline{\mathbf{u}}\end{array}\right),$$ E19

where ${\lambda }_p$ is the Lagrangian multiplier and both $R$ and $\tilde{R}$ are $1\times N_e$ matrices:

$$R=\left(\sqrt{6},0,\cdots ,0\right),\tilde{R}=\left(\frac{1}{N_e},\cdots ,\frac{1}{N_e}\right).$$ E20

The coefficient matrix on the left-hand side of Eq. (19) is the so-called reconstruction matrix [19, 20].

2.4 The spatial discrete formulation

We now derive the semi-discrete finite volume scheme based on Eqs. (2) and (8). Integrating over each subelement $T^{\left(m\left(k\right)\right)}$ on both sides of Eq. (2), we have

$${\int }_{T^{\left(m\left(k\right)\right)}}\frac{\partial \mathbf{u}}{\partial t}\mathit{dV}+{\int }_{T^{\left(m\left(k\right)\right)}}\left(A\frac{\partial \mathbf{u}}{\partial x}+B\frac{\partial \mathbf{u}}{\partial y}+C\frac{\partial \mathbf{u}}{\partial z}\right)\mathit{dV}=0,k=1,\cdots ,N_e.$$ E21

Using Eq. (8) and integration by parts yield

$${\int }_{T^{\left(m\left(k\right)\right)}}\frac{\partial \mathbf{u}}{\partial t}\mathit{dV}+{\int }_{\partial T^{\left(m\left(k\right)\right)}}{\mathbf{F}}^h\mathit{dS}=0,$$ E22

where $\mathit{dS}$ denotes the infinitesimal element in the face integral and ${\mathbf{F}}^h$ is the numerical flux, and we adopt the widely used Godunov flux [15, 19, 20, 23]

$$\begin{array}{rr}{\mathbf{F}}^h& =\frac{1}{2}T\left(A^{\left(m\left(k\right)\right)}+|A^{\left(m\left(k\right)\right)}|\right)T^{-1}\sum _{l=1}^{N_p}{\widehat{\mathbf{w}}}_l^{\left(m\right)}{\varphi }_l^{\left(m\right)}\\ & +\frac{1}{2}T\left(A^{\left(m\left(k\right)\right)}-|A^{\left(m\left(k\right)\right)}|\right)T^{-1}\sum _{l=1}^{N_p}{\widehat{\mathbf{w}}}_l^{\left(m_j\right)}{\varphi }_l^{\left(m_j\right)},\end{array}$$ E23

where $m_j$ is the index number of coarse tetrahedral element neighboring subelement $T^{\left(m\left(k\right)\right)}$ . The notation $\mid A^{\left(m\left(k\right)\right)}\mid$ denotes applying the absolute value operator of the eigenvalues given in Eq. (3), i.e.,

$$\mid A^{\left(m\left(k\right)\right)}\mid =R\mid \mathrm{\Lambda }\mid R^{-1},\mid \mathrm{\Lambda }\mid =\mathit{diag}\left(|s_1|,\cdots ,|s_9|\right),$$ E24

where $R$ is the matrix and its columns are made up of the eigenvectors associated with eigenvalues in Eq. (3), i.e.,

$$R=\left(\begin{array}{lllllllll}\lambda +2\mu & 0& 0& 0& 0& 0& 0& 0& \lambda +2\mu \\ \lambda & 0& 0& 0& 1& 0& 0& 0& \lambda \\ \lambda & 0& 0& 0& 0& 1& 0& 0& \lambda \\ 0& \mu & 0& 0& 0& 0& 0& \mu & 0\\ 0& 0& 0& 1& 0& 0& 0& 0& 0\\ 0& 0& \mu & 0& 0& 0& \mu & 0& 0\\ v_p& 0& 0& 0& 0& 0& 0& 0& -v_p\\ 0& v_s& 0& 0& 0& 0& 0& -v_s& 0\\ 0& 0& v_s& 0& 0& 0& -v_s& 0& 0\end{array}\right).$$ E25

And $T$ is the rotation matrix given by

$$T=\left(\begin{array}{llllll}{n}_x^2& s_x^2& t_x^2& 2n_x{s}_x& 2s_x{t}_x& 2n_x{t}_x\\ n_y^2& s_y^2& t_y^2& 2n_y{s}_y& 2s_y{t}_y& 2n_y{t}_y\\ n_z^2& s_z^2& t_z^2& 2n_z{s}_z& 2s_z{t}_z& 2n_z{t}_z\\ n_y{n}_x& s_y{s}_x& t_y{t}_x& n_y{s}_x+n_x{s}_y& s_y{t}_x+s_x{t}_y& n_y{t}_x+n_x{t}_y\\ n_z{n}_y& s_z{s}_y& t_z{t}_y& n_z{s}_y+n_y{s}_z& s_z{t}_y+s_y{t}_z& n_z{t}_y+n_y{t}_z\\ n_z{n}_x& s_z{s}_x& t_z{t}_x& n_z{s}_x+n_x{s}_z& s_z{t}_x+s_x{t}_z& n_z{t}_x+n_x{t}_z\end{array}\right),$$ E26

where $\left(n_x,n_y,n_z\right)$ is the normal vector of the face and $\left(s_x,s_y,s_z\right)$ and $\left(t_x,t_y,t_z\right)$ are the two tangential vectors. $T^{-1}$ denotes the inverse of $T$ .

Inserting Eqs. (23) into (22) and rewriting the result into a splitting form of easy computation in the reference system $\mathbf{\xi }$ , we have

$$\frac{\partial }{\partial t}{\overline{\mathbf{u}}}^{\left(m\left(k\right)\right)}\mid T^{\left(m\left(k\right)\right)}\mid +\sum _{j=1}^4{\mathbf{F}}_j^h=0$$ E27

with

$${\mathbf{F}}_j^h=T^j{A}^{\left(m\left(k\right)\right)}{\left(T^j\right)}^{-1}\mid S_j\mid \sum _{l=1}^{N_p}{F}_l^{-,j}{\widehat{\mathbf{w}}}_l^{\left(m\right)},m=m_j,$$ E28

and

$$\begin{array}{rr}{\mathbf{F}}_j^h& =\frac{1}{2}{T}^j\left(A^{\left(m\left(k\right)\right)}+|A^{\left(m\left(k\right)\right)}|\right){\left(T^j\right)}^{-1}\mid S_j\mid \sum _{l=1}^{N_p}{F}_l^{-,j}{\widehat{\mathbf{w}}}^{\left(m\right)}\\ & +\frac{1}{2}{T}^j\left(A^{\left(m\left(k\right)\right)}-|A^{\left(m\left(k\right)\right)}|\right){\left(T^j\right)}^{-1}\mid S_j\mid \sum _{l=1}^{N_p}{F}_l^{+,i,p}{\widehat{\mathbf{w}}}_l^{\left(m_j\right)},m\ne m_j,\end{array}$$ E29

where $S_j$ is the area of the $j$ -th $\left(j=\mathrm{1,2,3,4}\right)$ face of subelement $T^{\left(m\left(k\right)\right)}$ . $F_l^{-,j}$ and $F_l^{+,i,p}$ are the left flux matrix and the right state flux matrix, respectively, which are given by

$$F_l^{-,j}={\int }_{\partial {\left(T_E\right)}_j}{\varphi }_l\left(\mathbf{\xi }\left({\tilde{T}}^{\left(m\left(j\right)\right)},{\tilde{\mathbf{\xi }}}^{\left(j\right)}\left(\chi ,\tau \right)\right)\right)\mathit{\text{dχdτ}},j=\mathrm{1,2,3,4},$$ E30

$$F_l^{+,i,p}={\int }_{\partial {\left(T_E\right)}_j}{\varphi }_l\left(\mathbf{\xi }\left({\tilde{T}}^{\left(m\left(i\right)\right)},{\tilde{\mathbf{\xi }}}^{\left(i\right)}\left({\tilde{\chi }}^{\left(p\right)},{\tilde{\tau }}^{\left(p\right)}\right)\right)\right)\mathit{\text{dχdτ}},i=\mathrm{1,2,3,4};p=\mathrm{1,2,3}.$$ E31

where $\chi$ and $\tau$ are the face parameters. The transformation of the face parameters $\chi$ and $\tau$ to the face parameters $\tilde{\chi }$ and $\tilde{\tau }$ in the neighbor tetrahedron depends on the orientation of the neighbor face with respect to the local face of the considered tetrahedron. And the mapping is given in Table 2. For a given tetrahedral mesh with the known indices $i$ and $p$ , there are only 4 of 12 possible matrices $F^{+,i,p}$ per element [15, 20]. Comparing with the traditional FV method, the method with the splitting form described above has much less computations of face integrations. Note that only our proposed FV method can be written as a splitting form. Theoretical analysis shows our method can save about half computational time under the condition of the same number of elements [27].

$p$	1	2	3
$\tilde{\chi }$	$\tau$	$1-\chi -\tau$	$\chi$
$\tilde{\tau }$	$\chi$	$\tau$	$1-\chi -\tau$

Table 2.

Transformation of the face parameters $\chi$ and $\tau$ to the face parameters $\tilde{\chi }$ and $\tilde{\tau }$ .

2.5 The time discretization

Equation (27) is in fact a semi-discrete ordinary differential equation (ODE) system. In order to solve it formally, we denote the spatial semi-discrete part in Eq. (27) by a linear operator $L$ . Then Eq. (27) can be written as a concise ODE form:

$$\frac{d\mathbf{u}}{\mathit{dt}}=L\left(\mathbf{u},t\right).$$ E32

Traditionally, the classic fourth-order explicit RK (ERK) method

$$\begin{array}{l}{\mathbf{k}}^{\left(1\right)}=L\left({\mathbf{u}}^n,t^n\right),\\ {\mathbf{k}}^{\left(2\right)}=L\left({\mathbf{u}}^n+\frac{1}{2}\mathrm{\Delta }t{\mathbf{k}}^{\left(1\right)},t^n+\frac{1}{2}\mathrm{\Delta }t\right),\\ {\mathbf{k}}^{\left(3\right)}=L\left({\mathbf{u}}^n+\frac{1}{2}\mathrm{\Delta }t{\mathbf{k}}^{\left(2\right)},t^n+\frac{1}{2}\mathrm{\Delta }t\right),\\ {\mathbf{k}}^{\left(4\right)}=L\left({\mathbf{u}}^n+\mathrm{\Delta }t{\mathbf{k}}^{\left(3\right)},t^n+\mathrm{\Delta }t\right),\\ {\mathbf{u}}^{n+1}={\mathbf{u}}^n+\frac{1}{6}\mathrm{\Delta }t\left({\mathbf{k}}^{\left(1\right)}+2{\mathbf{k}}^{\left(2\right)}+2{\mathbf{k}}^{\left(3\right)}+{\mathbf{k}}^{\left(4\right)}\right)\end{array}$$ E33

can be applied to advance $\mathbf{u}$ from ${\mathbf{u}}^n$ to ${\mathbf{u}}^{n+1}$ . Here $\mathrm{\Delta }t$ is the time step. Now we use the low-storage version of ERK (LSERK) to solve Eq. (32):

$$\begin{array}{l}{\mathbf{u}}^{\left(0\right)}={\mathbf{u}}^n,\\ \left\{\begin{array}{l}{\mathbf{k}}^{\left(i\right)}=a_i{\mathbf{k}}^{\left(i-1\right)}+\mathrm{\Delta }\mathit{tL}\left({\mathbf{p}}^{\left(i-1\right)},t^n+c_i\mathrm{\Delta }t\right),\\ {\mathbf{p}}^{\left(i\right)}={\mathbf{p}}^{\left(i-1\right)}+b_i{\mathbf{k}}^{\left(i\right)},i=1,\cdots ,5,\end{array}\right.\\ {\mathbf{u}}^{\left(n+1\right)}={\mathbf{p}}^{\left(5\right)}.\end{array}$$ E34

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].

$i$	$a_i$	$b_i$	$c_i$
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

Table 3.

Coefficients for the low-storage five-stage fourth-order ERK method.

As to the stability condition, it is controlled by the Courant-Friedrichs-Lewy (CFL) condition [15, 19];

$$\mathrm{\Delta }t\le \frac{1}{2N+1}\frac{h_{\min }}{v_p},$$ E35

where $v_p$ is the $P$ wave velocity and $h_{\min }$ is the minimum diameter of the circumcircles of tetrahedral elements. This condition is a necessary condition for discrete stability, and a bit more restrictive form is actually used in numerical computations.

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

$${\mathbf{F}}^h=\frac{1}{2}T\left(A^{\left(m\left(k\right)\right)}+|A^{\left(m\left(k\right)\right)}|\right)T^{-1}\sum _{l=1}^{N_p}{\widehat{\mathbf{w}}}_l^{\left(m\right)}{\varphi }_l^{\left(m\right)},$$ E36

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 $\left(x_i,y_i,z_i\right)$ for $i=\mathrm{1,2,3,4}$ be the coordinates of a physical element. The transformation from $\xi -\eta -\zeta$ system to $x-y-z$ system is defined by

$$\left\{\begin{array}{l}x=x_1+\left(x_2-x_1\right)\xi +\left(x_3-x_1\right)\eta +\left(x_4-x_1\right)\zeta ,\\ y=y_1+\left(y_2-y_1\right)\xi +\left(y_3-y_1\right)\eta +\left(y_4-y_1\right)\zeta ,\\ z=z_1+\left(z_2-z_1\right)\xi +\left(z_3-z_1\right)\eta +\left(z_4-z_1\right)\zeta ,\end{array}\right.$$ E37

then the transformation from $x-y-z$ system to $\xi -\eta -\zeta$ system can be solved for $\xi ,\eta$ and $\zeta$ from Eq. (37) by the Cramer ruler, i.e.,

$$\xi =\frac{\mid J_1\mid }{\mid J\mid },\eta =\frac{\mid J_2\mid }{\mid J\mid },\zeta =\frac{\mid J_3\mid }{\mid J\mid },$$ E38

where

$$\mid J_1\mid =\left|\begin{array}{lll}x-x_1& x_3-x_1& x_4-x_1\\ y-y_1& y_3-y_1& y_4-y_1\\ z-z_1& z_3-z_1& z_4-z_1\end{array}\right|,\mid J_2\mid =\left|\begin{array}{lll}{x}_2-x_1& x-x_1& x_4-x_1\\ y_2-y_1& y-y_1& y_4-y_1\\ z_2-z_1& z-z_1& z_4-z_1\end{array}\right|,$$ E39

$$\mid J_3\mid =\left|\begin{array}{lll}{x}_2-x_1& x_3-x_1& x-x_1\\ y_2-y_1& y_3-y_1& y-y_1\\ z_2-z_1& z_3-z_1& z-z_1\end{array}\right|,\mid J\mid =\left|\begin{array}{lll}{x}_2-x_1& x_3-x_1& x_4-x_1\\ y_2-y_1& y_3-y_1& y_4-y_1\\ z_2-z_1& z_3-z_1& z_4-z_1\end{array}\right|.$$ E40

Note that $J$ is the determinant of the Jacobian matrix of the transformation being equal to six times the volume of the tetrahedron element $T^{\left(m\right)}$ .

The coordinate transformation from the second reference coordinate $\tilde{\xi }-\tilde{\eta }-\tilde{\zeta }$ to $\xi -\eta -\zeta$ system is defined by

$$\left\{\begin{array}{l}\xi ={\xi }_1+\left({\xi }_2-{\xi }_1\right)\tilde{\xi }+\left({\xi }_3-{\xi }_1\right)\tilde{\eta }+\left({\xi }_4-{\xi }_1\right)\tilde{\zeta },\\ \eta ={\eta }_1+\left({\eta }_2-{\eta }_1\right)\tilde{\xi }+\left({\eta }_3-{\eta }_1\right)\tilde{\eta }+\left({\eta }_4-{\eta }_1\right)\tilde{\zeta },\\ \zeta ={\zeta }_1+\left({\zeta }_2-{\zeta }_1\right)\tilde{\xi }+\left({\zeta }_3-{\zeta }_1\right)\tilde{\eta }+\left({\zeta }_4-{\zeta }_1\right)\tilde{\zeta },\end{array}\right)$$ E41

then the transform from $\xi -\eta -\zeta$ system to $\tilde{\xi }-\tilde{\eta }-\tilde{\zeta }$ system can be solved for $\tilde{\xi }-\tilde{\eta }-\tilde{\zeta }$ from Eq. (41) by Cramer ruler similarly. Denote

$$\mid \tilde{J}\mid =\left|\begin{array}{lll}{\xi }_2-{\xi }_1& {\xi }_3-{\xi }_1& {\xi }_4-{\xi }_1\\ {\eta }_2-{\eta }_1& {\eta }_3-{\eta }_1& {\eta }_4-{\eta }_1\\ {\zeta }_2-{\zeta }_1& {\zeta }_3-{\zeta }_1& {\zeta }_4-{\zeta }_1\end{array}\right|,$$ E42

which is the determinant of the Jacobian matrix of the transformation being equal to six times the volume of the subelement ${\tilde{T}}^{\left(m\left(k\right)\right)}$ for $k=1,\cdots ,N_e$ . In Eqs. (41) and (42), $\left({\xi }_i,{\eta }_i,{\zeta }_i\right)$ for $i=\mathrm{1,2,3,4}$ , denote the vertex coordinates of ${\tilde{T}}^{\left(m\left(k\right)\right)}$ in $\xi -\eta -\zeta$ system.

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 two-layered model with the inclined interface shown in Figure 3a. The range of the model is $x\in \left[-1.6\mathit{km},1.6\mathit{km}\right]$ and $z\in \left[-1.6\mathit{km}\mathrm{,1.8}\mathit{km}\right]$ . The surface of the model is uneven to imitate the real topography. The $v_p$ and $v_s$ velocities are $3000 \mathrm{m}/\mathrm{s}$ and $2000 \mathrm{m}/\mathrm{s}$ in the upper layer and $2400 \mathrm{m}/\mathrm{s}$ and $1600 \mathrm{m}/\mathrm{s}$ in the lower layer, respectively. The densities $\rho$ are $2200 \mathrm{kg}/{\mathrm{m}}^3$ and $1800 \mathrm{kg}/{\mathrm{m}}^3$ in the upper and lower layer, respectively. Figure 3b is the coarser triangular meshes for this model. A coarser version of the mesh is shown here as the finest mesh in computations cannot be seen clearly. The triangular meshes can fit the curve topography very well. Note that none triangular element crosses the interface. In computations the $P_4$ polynomial reconstruction is applied. The computational domain is meshed by 113472 coarse elements. Each coarse element is subdivided into 25 subelements further. So there are 2,836,800 fine elements totally. The time step is $\mathrm{\Delta }t=5\times {10}^{-5}\mathrm{s}$ . The source is located at $\left(x,z\right)=\left(\mathrm{0,0.2}\mathit{km}\right)$ with time history

Figure 3.

A two-layered model with curved surface topography (a) and the triangular meshes (b).

$$f\left(t\right)=-2\alpha \left(t-t_0\right)e^{-\alpha {\left(t-t_0\right)}^2},t_0=0.08,\alpha ={\left(\pi f_0\right)}^2,$$ E43

where $f_0=20\mathrm{Hz}$ is the main frequency. In order to simulate point source excitation, a spatial local distribution function defined by

$$G\left(\mathbf{x}\right)=\left\{\begin{array}{ll}\mathit{exp}\left(-7{‖\mathbf{x}-{\mathbf{x}}_0‖}_2^2/r_0^2\right),& {‖\mathbf{x}-{\mathbf{x}}_0‖}_2^2⩽r_0^2,\\ 0,& {‖\mathbf{x}-{\mathbf{x}}_0‖}_2^2>r_0^2,\end{array}\right.$$ E44

is applied, where ${\mathbf{x}}_0=\left(x_0,y_0,z_0\right)$ are positions of the source center. The source is added to the $u$ component; that is to say, all source terms except $g_7$ in Eq. (1) are all zero. Figure 4 is the snapshots of $u$ and $v$ components at propagation time $0.25\mathrm{s}$ . Figure 5 is the snapshots of $u$ and $v$ components at propagation time $0.30\mathrm{s}$ . We can see the $P$ wave and $S$ wave propagate toward out of the model. The reflected and transmitted waves due to the tilted physical interface are also very clear. These are the expected physical phenomena of wave propagation in elastic media.

Figure 4.

Snapshots of $u$ component (a) and $v$ component (b) at propagation time $\mathit{0.25}s$ .

Figure 5.

Snapshots of $u$ component (a) and $v$ component (b) at propagation time $\mathit{0.30}s$ .

Example 2. The second example is a cuboid model. The physical size of the model is $\left(x,y,z\right)\in \left[\mathrm{0,2}\mathit{km}\right]\times \left[\mathrm{0,2}\mathit{km}\right]\times \left[\mathrm{0,1}\mathit{km}\right]$ . The model and its unstructured tetrahedral meshes are shown in Figure 6. There are totally 836,612 coarse tetrahedrons to mesh the model. A coarser mesh is shown as the actual mesh in computations is too fine to see clearly. Each coarse tetrahedron is subdivided into $N_e=27$ subelements as we adopt $P_3$ polynomial reconstruction. The parameters for $\lambda$ , $\mu$ , and $\rho$ are ${10}^9\mathrm{Pa}$ , ${10}^9\mathrm{Pa}$ , and $1000\mathrm{kg}/{\mathrm{m}}^3$ . The time step in computations is ${10}^{-4}\mathrm{s}$ . The source is located in the center of the model with time history given by

Figure 6.

A cubic model and its unstructured tetrahedral meshes.

$$f\left(t\right)=sin\left(40\mathit{\pi t}\right)e^{-100t^2}.$$ E45

It is applied to the $u$ component. The 3D snapshots of $u$ , $v$ , and $w$ components at propagation time $0.42\mathrm{s}$ are shown in Figure 7. From these figures, we can clearly see two types of waves, i.e., the compressive wave and the shear wave. The splitting PML in nonconvolutional form is adopted here [32], and the boundary reflections are absorbed obviously and effectively. The message passing interface (MPI) parallelization based on spatial domain decomposition is applied. The CPU time for extrapolation 1000 time steps is about $33,310\mathrm{s}$ with 128 processors each with 2.6 GHz main frequency.

Figure 7.

The 3D snapshots of $u$ component (a), $v$ component (b), and $w$ component (c) at propagation time $\mathit{0.42}\mathrm{s}$ in a cuboid model. The source is located in the center of the model.

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 $x\in \left[\mathrm{0,2.0}\mathit{km}\right]$ , $y\in \left[\mathrm{0,3.5}\mathit{km}\right]$ , and $z\in \left[\mathrm{0,1.1}\mathit{km}\right]$ . The corresponding 3D mesh is shown in Figure 8b. A coarser version of the mesh is given as the actual mesh in computations is too fine to see clearly in the figure. The model is meshed with 210,701 relative coarse tetrahedral elements. Each coarse tetrahedron is subdivided into $N_e=64$ subelements as we adopt $P_4$ polynomial reconstruction, and thus there are 13,484,864 fine elements totally. The time step $\mathrm{\Delta }t$ is ${10}^{-4}\mathrm{s}$ . The source is situated at $\left(x_0,y_0,z_0\right)=\left(750m,1300m,300m\right)$ with the same time history in Eq. (45). The media velocities of $v_p$ and $v_s$ are $v_p=3000\mathrm{m}/\mathrm{s}$ and $v_s=2000\mathrm{m}/\mathrm{s}$ . The MPI parallelization based on spatial domain decomposition is applied. The nonconvolutional splitting PML [32] is adopted. The 3D snapshots of $u$ , $v$ , and $w$ components at propagation time $0.80\mathrm{s}$ are shown in Figure 9. The CPU time for extrapolation 10,000 time steps is $100,449\mathrm{s}$ with 256 processors each with 2.6 GHz main frequency. From Figure 9, we can see clearly the propagation of $P$ wave and $S$ wave.

Figure 8.

A real 3D model with complex topography. (a) model and (b) unstructured tetrahedral meshes.

Figure 9.

3D snapshots of $u$ , $v$ , and $w$ components at propagation time $\mathit{0.80}\mathrm{s}$ in a real 3D model. The results are obtained by the method in this chapter with $P_4$ reconstruction. (a) u component, (b) v component, (c) w component.

4. Conclusions

A new efficient high-order 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 high-order polynomial reconstruction on the stencil. The resulting reconstruction matrix remains unchanged for all tetrahedrons and can be pre-computed 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 large-scale 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.