Abstract
This study concentrates on transient multiphysical wave problems for simulating seismic waves. The presented models cover the coupling between elastic wave equations in solid structures and acoustic wave equations in fluids. We focus especially on the accuracy and efficiency of the numerical solution based on higher-order discretizations. The spatial discretization is performed by the spectral element method. For time discretization we compare three different schemes. The efficiency of the higher-order time discretization schemes depends on several factors which we discuss by presenting numerical experiments with the fourth-order Runge-Kutta and the fourth-order Adams-Bashforth time-stepping. We generate a synthetic seismogram and demonstrate its function by a numerical simulation.
1. Introduction
Several formulations exist for modeling the seismic vibrations as interaction between acoustic and elastic waves. Typically, the displacement is solved in the elastic structure. The fluid can be modeled using finite element formulations based on fluid pressure, displacement, velocity potential, or displacement potential [1]. Two approaches, in which the displacement is solved in the elastic structure, predominate in modeling the interaction between acoustic and elastic waves. Expressing the acoustic wave equation by the pressure in the fluid domain leads to a nonsymmetric formulation (see, e.g., [2, 3]), while using the velocity potential results in a symmetric system of equations (see, e.g., [4–7]). Recently, a velocity-strain formulation has been considered by Wilcox et al. [8]. Essentially, the difference between the different models is only in the choice of the variables presenting the wave propagation in different domains. Nevertheless, from the numerical point of view, this choice determines the features of the coupled problem. Hence, it also gives the guideline for using appropriate discretization and solution methods.
Discretization methods play a crucial role in the efficiency. The key factor in developing efficient solution methods is the use of high-order approximations without computationally demanding matrix inversions. We attempt to meet these requirements by using the spectral element [9] method (SEM) for space discretization. The SEM was pioneered in the mid 1980s by Patera [10] and Maday and Patera [11], and it combines the geometric flexibility of finite elements with the high accuracy of spectral methods. The method is widely used for simulating seismic waves, in both frequency and time domains (see, e.g., [12, 13]). We have applied it to time-harmonic acoustoelastic equations in [14, 15], while this paper concentrates on the accuracy and efficiency of the time domain solutions.
In time domain simulations, the efficiency of the method depends also strongly on the time discretization. The second-order central finite difference (CD), or leap-frog, scheme is a low-order method which is easy to implement and thus in general use for discretizing in time domain. This method is used with the SEM by Komatitsch et al. [16], and a modification for utilizing larger time steps in the fluid than in the solid media is presented by Madec et al. [17]. Since the method is only of second-order accuracy, higher-order time discretizations are needed for efficient computer simulations based on higher-order space discretization. The fourth-order central finite differences are not usable at this stage, since the scheme is unconditionally unstable. Instead, a second-order accurate time discretization scheme can be converted to a fourth-order accurate one, as is done by Tarnow and Simo [18]. This symplectic scheme is compared with the CD method by Nissen-Meyer et al. [19]. The importance of higher-order time discretizations, in the field of seismology, has recently been considered by, for example, De Basabe and Sen [20], Peter et al. [21], and Liu et al. [22] and references therein.
Lax and Wendroff [23] introduced a higher-order generalization of the finite difference scheme. The method has been applied for acoustic wave equations to provide the fourth-order accuracy, with the finite difference space discretization by Dablain [24] and with the spectral element space discretization by Cohen and Joly [25]. However, this approach does not fit to a case, in which an absorbing boundary condition is used for truncating the computational domain, unless velocity-stress or velocity-displacement formulation is used [26].
By Kubatko et al. [27], Runge-Kutta time discretization methods were used in conjunction with discontinuous Galerkin (DG) finite element spatial discretization, and Antonietti et al. [28] have applied the scheme for elastic wave propagation problems. The second- and third-order Runge-Kutta methods are presented, with a fourth-order Hamiltonian-based space discretization for acoustic and elastic waves in separate domains by Ma et al. [29] and with -adaptive discontinuous Galerkin method for simulating tsunami propagation by Blaise and St-Cyr [30], respectively. Since there were no research results considering the accuracy and efficiency of the numerical simulation of fluid and solid waves, in which the fourth-order Runge-Kutta (RK) time discretization is combined with higher-order space discretization, we have applied it for acoustic [31] and elastic [32] waves. The method is used with the eighth-order finite difference space discretization for solving an elastodynamic problem also by Martin et al. [33].
With respect to the time step , the CD method is second-order accurate, while the RK method is fourth-order accurate. Although the computational effort of the RK method is approximately four times that of the CD scheme at each time step, the results of using higher-order time-stepping scheme were promising [31, 32]. Both methods lead to an explicit time-stepping scheme. The drawback is that the schemes need to satisfy the stability condition, which limits the length of the time step. For acoustic waves, Yang et al. [34] show that a fourth-order Runge-Kutta method is more accurate for time discretization than the fourth-order Lax-Wendroff method.
In this paper, we apply the CD and the RK methods to multiphysical wave problems. We also present the fourth-order Adams-Bashforth (AB) scheme to discuss the efficiency factors of the higher-order time discretization schemes. In what follows, we first consider the mathematical models and boundary conditions in Section 2. Further, we discretize the models with respect to space and time in Sections 2.2 and 2.3. Numerical examples are considered for the validation of the accuracy of the approximations and for demonstrating the synthetic seismogram in Section 3. The concluding remarks are presented in Section 4.
2. Numerical Models
In what follows, we concentrate essentially on two multiphysical linear models for simulating the transient propagation of seismic waves. In both models, the domain is divided into the solid part and the fluid part (see Figure 1). First, we present the model for elastic deformations in and use it as a starting point for deriving the different models presenting the propagation of acoustic vibrations in the fluid domain . Then, we apply the coupling conditions to the models and discretize the problems. We see that, despite the fact the difference between the models is only in the choice of the variable presenting the wave propagation, it affects the structure of the coupled problem at the discrete level. For the sake of the efficiency of the solving process and the accuracy of the numerical solution, the numerical methods, relying on the model-specific features, are presented.
2.1. Coupled Problems
The deformation of the elastic structure is modeled by the displacement in the solid domain , such thatwhere is the density of the structure and is the source function. We concentrate on isotropic materials, for which the number of essential material constants reduces to the two Lamé parameters, and , where is the Young modulus describing the stiffness of the solid and is the Poisson ratio presenting the compressibility of the solid as the ratio of lateral to longitudinal strain in a uniaxial tensile stress, such that . Hence, the stress tensor can be presented as , where is the identity matrix. In general, we assume the medium to be heterogeneous. Thus, the partial derivatives in apply to and as well as to the displacement. The speed of pressure waves (P-waves) and shear waves (S-waves) are presented as functions of the Lamé parameters and density. The P-waves move in a compressional motion, while the motion of the S-waves is perpendicular to the direction of wave propagation [35].
In principle, the wave equation in fluid media can be derived as a special case of (1), by taking into account the fact that there are no S-waves in fluids. That is because fluids are inviscid and, thus, have no internal friction and can not support shear stresses. Thus, taking the divergence of (1) leads to modeling acoustic waves in fluid domain by the pressure field (see, e.g., [16, 36]),where is the density of fluid, is the speed of the wave, depends on the vector-valued displacement in the fluid domain, and is the source function. By using the velocity potential , the corresponding model iswhere is the source function (see, e.g., [4, 5, 17]).
The equations need to be completed by the initial and boundary conditions to get a well-posed and physically meaningful problem. In this paper, we use the first-order absorbing boundary condition [37] for truncating the exterior domain on the boundaries and . On the interface between fluid and solid domains, the normal components of displacements and forces are balanced. That is presented, depending on the choice of variable in the fluid domain, byorwhere and are the outward pointing unit normal vector on the boundary , with respect to the solid and fluid domain, respectively. The interface conditions coupling the two domains are needed to model the following fact: when a wave coming from the fluid domain confronts the elastic domain, it is not totally reflected, but part of it passes to the elastic domain and turns to elastic vibrations. The analogous action is seen when the elastic wave propagates to the fluid, although usually the reflections from fluid to solid are minor when compared with the reflections from solid to fluid. To be more precise, the magnitude of the reflection depends on the difference between the densities and the wave speeds of the materials. Thus, the reflections are more significant from comparatively stiff structures than for more flexible obstacles.
For the weak formulation of the system consisting of (1)-(2) and (4) we introduce the function spaces and . By multiplying (1) with any test function in the space and (2) with any test function in the space , using Green’s formula, and substituting the boundary conditions, we get the following weak formulation: find satisfying for any andwhere and are the source term acting on the artificial absorbing boundaries. Respectively, we can define the weak formulation for the system of (1), (3), and (5).
2.2. Spatial Discretization
The spectral element method is obtained from the weak formulation of the coupled problem by restricting the problem presented in the infinite-dimensional spaces into finite-dimensional subspaces. Thus, in order to produce an approximate solution for the problem, the given domain is discretized into a collection of quadrilateral elements , , such that . Each element is associated with four nodes. For the discrete formulation, we define the reference element and invertible affine mappings such that . The boundary of is a union of its four edges, which are the images of the four edges of the reference element, obtained by the mapping . Respectively, the vertex nodes of are the images of the four vertices of the reference element, obtained by the mapping . These properties are essential for the mappings to be sufficiently smooth. Each of elements is individually mapped to the reference element, and we make use of the affine mapping to make transformations from the physical domain to the reference domain and vice versa. The mapping between the reference element and the th element is defined such that , where and are the Gauss-Lobatto points in the reference element (see, e.g., [9]).
After the spatial discretization by the spectral elements, the semidiscrete form of the coupled problem is stated aswhere is the global block vector containing the values of the variables in both fluid and solid domains at time at the Gauss-Lobatto points of the quadrilateral mesh. The variable which we use in the solid domain is the displacement . Thus, the number of degrees of freedom (DOF) in the solid domain, , is in the two-dimensional domain twice the number of discretization points in the solid domain . The total number of degrees of freedom, that is, , depends on the variable of the fluid domain. If the fluid domain is modeled by using pressure or velocity potential , we have scalar values at each spatial discretization point, and the number of degrees of freedom in the fluid domain, expressed as , is equal to the number of discretization points in the fluid domain .
In practice, the discrete counterparts of the variables are approximated as linear combinations of the corresponding nodal values and the spectral element basis functions , , and , , which are higher-order Lagrange interpolation polynomials [9]. In the case of the formulation with pressure and displacement, entries of the matrices , , and and the right-hand side vector are given by the formulaswhere the matrix blocks and the -dimensional right-hand side vector corresponding to the fluid domain arewhere . Respectively, the block matrices and the -dimensional vector representing the elastic waves are which have the components where . The matrices arising from the coupling between acoustic and elastic wave equations are and , for which it holds thatFor , and , whereas, for , and .
The computation of the elementwise matrices and vectors involves the integration over the elementwise subregions. Evaluating these integrals analytically is usually complicated or even impossible. That is why a numerical integration procedure is used. In practice, we replace the integrals by finite sums, in which we use Gauss-Lobatto weights and nodal points. The values of these sums are computed element by element with the Gauss-Lobatto integration rule. Collocation points are now the nodes of the spectral element. All but one of the shape functions will be zero at a particular node. Thus, for , and implying that the matrices and are diagonal. Furthermore, the matrix is a lower triangular block matrix with diagonal blocks. Thus, the inverse of the matrix is also a lower triangular block matrix with diagonal blocks, and explicit time-stepping with central finite differences requires only matrix-vector multiplications.
In practice, the stiffness matrix is assembled once at the beginning of the simulation. It is stored by using the compressed column storage including only the nonzero matrix elements. The other options would have been using a mixed spectral element formulation [38].
In the case of the formulation with velocity potential and displacement, the entries of the matrices , , and and the right-hand side vector are given by the formulas where the matrix blocks and the -dimensional right-hand side vector corresponding to the fluid domain arefor . The block matrices and the -dimensional vector representing the elastic waves are exactly the same as in the nonsymmetric formulation. The matrices arising from the coupling between acoustic and elastic wave equations are and , for which it holds thatFor , and , whereas, for , and .
Using the vector-valued displacement in the fluid domain doubles the number of degrees of freedom in the fluid domain. In other words, . Therefore also the memory consumption for computing the pure displacement-displacement interaction is higher than in the case of the other couplings considered. Furthermore, if the formulation with displacement in both domains is considered, the divergence of the variable in the fluid domain is involved in the weak formulation, and we would need a scheme that approximates the functions in better than the spectral element method does. For instance, Raviart-Thomas finite elements, which are used for acoustic wave equations [39], could be used for this purpose. Furthermore, the Raviart-Thomas elements are not a good choice for discretizing the solid domain. Nevertheless, the coupling conditions should be satisfied at the interface of the two domains. That is, we would need to change the discretization approach in the solid domain as well, to make the degrees of freedom coincide or fulfill the coupling conditions, for example, by using Lagrange multipliers [40].
2.3. Time Discretization
After dividing the time interval into time steps, each of size , applying the appropriate time discretization into semidiscretized form (7), and taking into account the initial conditions, we obtain the matrix form of the fully discrete state equation. Despite the use of higher-order space discretization, the time-stepping for transient wave equations is usually performed by low-order methods like the central finite difference (CD) scheme. By proceeding this way, in the case of the nonsymmetric formulation, at each time step we compute first the displacement and then the pressure from the equationswith the initial conditionswhere , , , and are the vectors , , , and at . Because the matrix sums and are diagonal, their inverses are obtained simply by inverting each diagonal element. Respectively, the state equation for the formulation with velocity potential and displacement iswhere , , , and are the vectors , , , and at . In this case, and need to be solved simultaneously, and we need to invert at each time step. If the boundaries of the computational domain consist only of horizontal and vertical lines, the coefficient matrix needed for inversion at each time step can be implemented as a block matrix consisting of diagonal blocks. Then, and can be solved simply by using matrix-vector multiplications from the formulas where and are diagonal matrices and
Remark 1. It would also be possible to uncenter one of the two first-order derivatives in time to uncouple the problem. That approach would involve first-order difference approximations and introduce some dissipation. That is why we neglect deeper considerations of the uncentered schemes.
Sate equation (7) can be presented as a system of differential equations , where is a vector of time-stepping variables and and the function . To this modified form, we can apply the fourth-order Runge-Kutta method, which is a Taylor series method. In general, the Taylor series methods keep the errors small, but there is the disadvantage of requiring the evaluation of higher derivatives of the function . The advantage of the Runge-Kutta method is that explicit evaluations of the derivatives of the function are not required, but linear combinations of the values of are used to approximate . In the fourth-order Runge-Kutta method, the approximate at the th time step is defined aswhere contains the global block vector , including the values of the variables in both the fluid and the solid domain at the th time step, and its derivative at time , . The initial condition is given by , and , , are the differential estimates as follows:In other words, in order to get differential estimates (23), the function is evaluated at each time step four times, and then the successive approximation of is calculated by formula (22).
If the matrix is diagonal, as is in the formulation with the velocity potential, the only matrix inversion needed in time-stepping is computed simply by inverting each diagonal element in the matrix . This requires only floating point operations, which is the number of diagonal elements in the matrix and known as the number of degrees of freedom in the space discretization. Since the matrix contains only diagonal blocks and coupling terms, the operation count of the matrix-vector product is of order . In the matrix-vector multiplication involving the sparse stiffness matrix , only nonzero matrix entries are multiplied, which requires the order of operations, where is the order of the polynomials in the spectral element basis. Besides these, additions and multiplications are needed for a single evaluation of the function . According to (22), the computation of needs floating point operations. Thus, the computational cost for each time step of the state equation is of order also with the RK time-stepping. Although the computational cost is of the same order for both the CD and the RK time-steppings, the number of floating point operators needed for the RK is nearly four times that of the CD.
Next, we make an effort to decrease the computing time and still maintain the high accuracy provided by higher-order time discretizations. For this purpose, we present the fourth-order Adams-Bashforth (AB) method based on approximating the functions by interpolating polynomials. It gives the solution at the th time step aswhere contains the vector and its derivative at time , . Hence, it is a multistep method that requires information at four previous time steps implying that another method is needed for starting the time marching. At this stage, we utilize the fourth-order Runge-Kutta method for computing the start-up values , , , and to initialize the multistep method.
3. Numerical Examples
In this section, the computational accuracy and efficiency of solving seismic wave problems are considered. The focus of Section 3.1 is on different formulations and implementation processes. In Section 3.2, the efficiency is further analysed by using higher-order time discretization. Finally, simulated seismograms and two-dimensional illustrations presenting the seismic wave motion are given in Section 3.3. The numerical experiments presented in this section are carried out on an AMD Opteron 885 processor at 2.6 GHz.
3.1. Comparison between Symmetric and Nonsymmetric Formulations
We illustrate the computational cost of both symmetric and nonsymmetric formulations by solving a time-dependent problem both expressing the acoustic wave equation by the pressure and by using the velocity potential in the fluid domain. The right-hand sides and initial conditions are defined to satisfy the analytical solution , , and .
The problem is solved in a domain, which consists of the solid part and the fluid part (see Figure 1). We use square-element meshes with mesh step size , and the element order is increased in both parts of the domain from 1 to 5. The meshes are matching on the coupling interface set at for . On the other boundaries we have the absorbing boundary conditions. The material parameters in the fluid domain are and . In the solid domain, we use the values , , and . The angular frequency is the same for both media, and we set the propagation direction by the vector . The time interval is divided into 400 steps, each of size , to guarantee the stability condition also for the higher element orders.
When the pressure formulation is considered, the time marching involves only matrix-vector multiplications, and actual matrix inversions are not needed. In the implementation of the velocity potential formulation, a “one-shot” method is required for solving the linear system including and at each time step . In that case, the matrix which is needed to be inverted is stored either as a band matrix or by using the compressed column storage including only the nonzero matrix elements.
The Lapack LU decomposition routines dgbtrf and dgbtrs use the band matrix storage mode. With these routines the memory and CPU time consumption increases rapidly when the element order is increased (see Figure 2). The SuperLU library routines perform an LU decomposition with partial pivoting, and the triangular system solves through forward and backward substitution. At this stage, we also utilize the sparsity of the matrix by using the compressed column storage mode. Since all the nonzero elements are not near the diagonal of the matrix, the compressed column storage mode requires less storage and computing operations than the band storage mode. That is why the SuperLU library gives a less demanding procedure for solving the linear system than the Lapack library. This is because the matrices arising from the space discretization and including the coupling terms have, in general, a sufficiently large bandwidth. Thus, a remarkably larger amount of memory is needed for storing the sparse coefficient matrix in the band matrix form used in conjunction with Lapack routines than in the compressed column storage utilized with the SuperLU. Consequently, less time is needed when fewer elements are employed in the solution procedure with the SuperLU. The performance of the Lapack routines could be improved, for instance, by using an appropriate node numbering of the mesh.
(a) CPU time with respect to the element order
(b) Memory consumption with respect to the element order
We conclude that the computational efforts are of the same order of magnitude whether the problem in the fluid domain is solved with respect to pressure or whether we use velocity potential formulation in conjunction with the linear solver provided by the SuperLU library.
The comparison between numerical and analytical solution shows that in both media the accuracy improves when the element order grows until a certain error level is reached (see Figure 3). This error level, shown as a horizontal line in Figure 3, reflects the error level of time discretization. Since we use a fixed number of time steps at all element orders, the error of time discretization becomes dominant for higher-order elements. Hence, finer time steps or higher-order time discretizations are needed in conjunction with higher-order elements.
(a) Solid domain
(b) Fluid domain
In principle, both symmetric and nonsymmetric formulation should lead to the same order of accuracy at each element order. The results obtained in the solid domain and depicted in Figure 3(a) are perfectly in balance with this hypothesis. However, in Figure 3(b) we see that in the fluid domain higher accuracy is obtained with the velocity potential than with the pressure formulation at each element order. In the case of pressure formulation, the interface condition is derived by differentiating the displacement component twice with respect to time. From the physical point of view, some information is lost in that procedure, and linear growth of the displacement with respect to time is not eliminated at the interface. However, the solutions of the pressure formulation converge towards the analytical solution. Although the pressure formulation gives less accurate results than the velocity potential formulation, it is still a plausible choice at some point. For instance, the implementation process can be hastened when the preexisting solvers of acoustic and elastic problems can be harnessed in the implementation, and no additional linear solvers are needed.
3.2. Higher-Order Time Discretizations
We demonstrate how the efficiency of the method can be improved by using the fourth-order Runge-Kutta scheme instead of the central finite difference time discretization. In principle, the central finite difference scheme is second-order accurate, while the fourth-order Runge-Kutta method is fourth-order accurate. However, the error of space discretization limits the accuracy of the overall time-stepping scheme as well, and the stability issue restricts choosing a feasible length of the time step. To better illustrate the error of time and space discretization, we continue with the example presented in Section 3.1, except that we change the mesh step size to and use various time step lengths to observe stability and accuracy issues. The number of time steps needed for stability is first determined numerically by using time steps per time period, for , until a stable solution is achieved. From these results, we can define stability constant for each element order such thatThe stability conditions corresponding to the largest stable time step are given in Table 1. The stability region seems to be exactly the same with both the CD and the RK time-stepping for the element orders and differs only slightly for the higher-order elements. The results reflect the well-known CFL condition and are in a good agreement with the experiments presented for acoustic waves with the CD time discretization by Cohen [9].
We start the computations with the largest stable time step and then repeatedly add the number of time steps by 300, until their number is larger than 3000. Proceeding in this way, for each element order we achieve a series of numerical results with various lengths of the time step. The maximum errors between the numerical and the analytical solution with respect to are computed as -norms. Accuracy of the numerical solution is shown in Figure 4 with respect to the ratio between the time step and the mesh step size for both the CD and the RK time-steppings with five element orders . Every curve represents computations with a particular polynomial order which has a characteristic discretization error. Naturally, the order of the space discretization error decreases when higher-order elements are used. It is worth mentioning that errors are somewhat smaller in the solid domain than in the fluid domain.
(a) Solid domain
(b) Fluid domain
It is seen that with the element orders and both time-stepping schemes give the same accuracy even when sufficiently large time steps are used. Moreover, this is the accuracy of the space discretization since the error of spatial discretization dominates with low-order elements. In other words, the spatial error is larger than the temporal error, and the maximum error with respect to the length of the time step is not decreasing significantly even if smaller time steps are used.
Naturally, for each element order the solution with the RK time-stepping is at least as accurate as the one computed with the CD time-stepping. When higher-order elements are used, results computed with the RK time discretization are more accurate than the ones computed with the CD time discretization. Depending on the accuracy of the time discretization, the error of temporal discretization might be dominating with large time steps. This is shown especially in the case of the CD time discretization with space discretization orders . In principle, also with these higher-order elements, the error curves turn to vertical lines, reflecting the accuracy of the space discretization, when the length of the time step is refined enough. In the case of the RK time discretization, the error of space discretization is dominant even when long time steps are used. With the RK time discretization, very fine time steps are needed only for in the solid domain to achieve the error level of spatial discretization.
With the help of the results depicted in Figure 4 we can extrapolate that very fine time steps are needed with the CD time discretization to keep the temporal error smaller than the spatial error when higher-order elements are involved. For , we would already need thousands of time steps with the CD time discretization to get the same accuracy as with the RK time discretization with 350 time steps. Tens of thousands of time steps are required to conceal the temporal error with the CD time discretization with . Respectively, for we need hundreds of thousands of time steps to achieve the error level of space discretization. The values of , satisfyingand concealing the temporal error, are reported for different element orders in Table 2.
Since refining the time steps means more evaluations, we are also interested in measuring the CPU time needed in these simulations. The numerical results about CPU time consumption are seen in Figure 5. To conceal the temporal error for and , less computation time is needed with the CD time discretization than with the RK time discretization. With higher-order elements, more remarkable time saving occurs by using the RK time discretization to attain the accuracy of the space discretization. Although the difference in time consumption is not significant with a sufficiently small number of degrees of freedom, it can play an important role in large-scale real-life applications.
(a) Solid domain
(b) Fluid domain
Next, we consider the Adams-Bashforth (AB) time-stepping scheme, with which only one evaluation of the function is needed at each time step. The values are stored to be used in the four following time steps. Hence, the memory consumption is of the same order of magnitude as when the fourth-order RK method is utilized. In theory, employing the fourth-order AB time-stepping scheme, instead of the fourth-order RK method, saves time. Nevertheless, the stability region is smaller for the fourth-order AB than the fourth-order RK method. The largest stable time steps with the AB method, determined numerically in the same way as is done in the previous example, are reported in Table 3. This practical realization shows that the length of the time step that guarantees the stability conditions is much smaller with the AB method than with the RK method. The computations are continued by carrying out the simulations with the AB time discretization with smaller time steps. That is, for the element order the addition of time steps is repeated until the number of time steps is larger than 3000. The results are presented in Figure 6 as accuracy with respect to CPU time consumption. For comparison, the results of the RK time discretization from the previous example are also depicted there. Clearly, stability restrictions deteriorate the efficiency of the AB time-stepping scheme. The CPU time used for the computation is in favor of the RK time discretization; time that is more than an order of magnitude longer is consumed to get the same accuracy with the AB time discretization than with the RK time discretization. Thus, the efficiency of the method can not be improved by utilizing the AB time discretization. Furthermore, the AB time discretization requires the solutions of the earlier time steps to be stored and, hence, is not memory efficient.
(a) Solid domain
(b) Fluid domain
3.3. Simulated Seismogram
In the last example, we consider a Runge-Kutta time-simulation for the velocity potential formulation with the element order . We use the same computational domain as in the previous examples with matching square-element meshes with step size . The material parameters are , , , , and . A point source is set on the fluid domain, such thatIt generates an impulse and is typical for modeling seismic waves related to earthquakes or explosions (see, e.g., [29, 41, 42]). The other source functions, , , and, , as well as the initial conditions, , , and , are zero-valued. The simulation has been run from to with 3100 time steps, each of size . The results, observed as the values of the velocity potential at the point , are presented in Figure 7. Respectively, the horizontal and vertical components of the displacement at the point are seen in Figure 8. That is, we have simulated results corresponding to measurements obtained by a receiver, such as a seismometer, placed at a particular place. Further, we present in Figure 9 a two-dimensional snapshot of both the domains at time step 1450.
(a) Horizontal component
(b) Vertical component
(a) Displacement amplitude in
(b) Velocity potential in
4. Conclusions
We considered the spectral element space discretization for time-dependent equations considering acoustic and elastic wave propagation and their interaction. It is possible to construct the spectral element formulation such that it results in a global mass matrix that is diagonal by construction, which leads to efficient implementation. This is an advantage compared with the classical finite element method.
The temporal discretization for time-dependent equations was made by the second-order accurate central finite difference scheme or the fourth-order accurate Runge-Kutta or Adams-Bashforth approaches. The performances of the different time discretization schemes were compared numerically. We found out that the Adams-Bashforth scheme has such strict stability conditions that it is not a feasible choice at this stage. The CPU time used for the computation is in favor of the RK time discretization; namely, time that is more than one order of magnitude longer is consumed to get the same accuracy with the AB time discretization than with the RK time discretization.
We can conclude that, to make good use of higher-order elements, also the time discretization should be done with a higher-order scheme. As a rule of thumb we can say that the efficiency of the overall method suffers from the error of time discretization if the order of the element is greater than the order of the time discretization method used. The second-order central finite difference time discretization method is efficient with finite elements, but when high accuracy is needed, it is best to use efficient higher-order time discretization methods such that the Runge-Kutta scheme.
Competing Interests
The author declares that there are no competing interests regarding the publication of this paper.
Acknowledgments
The author would like to thank Professor Jari Toivanen and Professor Tuomo Rossi for useful discussions. The research was supported by the Academy of Finland, Grants 252549 and 259925.