Frequency-domain Elastic wave Simulation Based on the Nonoverlapping Domain Decomposition Method

A new wave simulation technique for the elastic wave equation in the frequency domain based on a no overlapping domain decomposition algorithm is investigated. The boundary conditions and the finite difference discrimination of the elastic wave equation are derived. The algorithm of no overlapping domain decomposition method is given. The method solves the elastic wave equation by iteratively solving sub problems defined on smaller sub domains. Numerical computations both for homogeneous and inhomogeneous media show the effectiveness of the proposed method. This method can be used in the full-waveform inversion.


Introduction
The numerical solution of wave equations is a very important problem.The applications can be found in electromagnetic, acoustic wave propagation, seismic inversion.The wave equations can be solved in either the time or the frequency domain.When applying the Fourier transformation with respect to time to the acoustic wave equation, the Helmholtz equation is obtained.The Helmholtz equation has important applications in underwater and seismic imaging [1,2].In the geophysical frequency-domain inversion one needs to do forward modeling which means solving the Helmholtz equation.During the inversion process, the synthetic model is continuously updated until a convergence is reached.Thus finding an efficient numerical method to solve the Helmholtz equation is important.For a finite element or finite difference discrimination of the 2D Helmholtz, it is necessary to use large number of grid points to handle large domains and high wave numbers.Moreover, the requirement of ten points per wavelength at least leads to very large complex linear systems.In seismic fullwaveform inversion, the number of grid points grows very rapidly in order to maintain good accuracy.The error is proportional to where k is the wave number, p is the order of discrimination and h is the grid size [3,4].The linear system becomes extremely large and highly indefinite.This makes the problem even harder to solve.Direct methods easily suffer from inacceptable computational work.The iterative methods for the Helmholtz equation have been an active research field since the 1980s.Since then many attempts have been spent to develop a powerful iterative method for solving it, see, for instance, [2,[5][6].
The domain decomposition method (DDM) is an effective technique for solving large-scale problems [7][8][9][10].The idea is to split the computational domain  into several smaller m sub domains i , and solve a sequence of similar sub problems on these subdomains.The boundary conditions are adjusted iteratively the transmission conditions between adjacent subdomains.The number and size of sub domains can now be chosen so as to enable direct methods to solve the sub problems.Several classes of DDM exist, e.g.multiplicative and additive Schwartz methods [7][8][9][10].
Absorbing boundary conditions (ABCs) are required in numerical simulations of unbounded problems.They are enforced at edges of a computational domain to absorb outgoing waves and thereby model an unbounded region.A number of ABCs have been introduced for use in finite-difference simulation of elastic wave propagation.There exist several ABC methods, for instance, the sponge layer method [11], the paraxial approximation method [12] and the perfectly matched layer method [13].The approach based on a paraxial approximation of the elastic wave equations was introduced by Clayton and Engquist [12].It is often adopted in frequency-domain simulation.The most recent advance in ABCs theory is the so-called perfectly matched layer (PML) technique proposed by Berenger [13] which offers an almost complete elimination of boundary reflections.The PML method has been widely used in wave simulation.
The DDM solving procedure for the Helmholtz has been investigated; for example, see [14][15][16].In this paper, we focus on the elastic wave simulation in the frequency domain rather than the acoustic wave based on the no overlapping domain decomposition method.Here we use the no overlapping DDM algorithm to solve the elastic wave equation in the frequency domain.

Finite Difference Discretization
The 2-D, second-order, frequency-domain elastic wave equation in a homogeneous, isotropic, and source-free media consists of the two coupled equations where 2 f    is the angular frequency,  is the density, and  and  are Lame parameters.The wavefield variables and are, respectively, the horizontal and vertical components of the Fourier transformed displacements.In a more general formulation, one also included a source term in the right-hand side of (1) or (2).

u v
We wish to find the solution to system (1)-( 2) for a finite model domain, , but without reflected energy from the domain boundary.For simplicity we use a simple absorbing boundary condition firstly proposed by Clayton and Engquist [11 where n is the direction normal to the absorbing model boundary, and are the velocities of compressive wave and shear wave, respectively.The condition ( 3) is only perfectly absorbing when the out-going wave impacts the boundary at normal incidence.We discredited (1)-( 3) using an implicit finite-difference scheme based on the 2nd order accurate centered-difference operators.We assume the computational domain is rectangular.A 2D regular mesh with constant grid spacing is used for deriving computational schemes.The difference schemes for (1)-( 2) are the following: ) 4 The difference schemes for the absorbing boundary condition (3) can be written as ,0 , , , , where x  and z  are the spatial steps in x and z directions respectively.The difference schemes ( 7)-( 10) are used at the top, bottom, left and right boundary respectively.Combining all the difference schemes above, we have the following concise form where A is the known matrix, ij A can be considered the sub matrixes of A , 1 and 2 b are the known vectors.They can be obtained based on the difference schemes above.For example, the elements of b 11 A can be determined by the following expressions: ( ( 1) , ( 1) ) ( ( 1) , ) , ( ( 1) , ( , ; r is the grid ratio / r x   z ; x and z are the discretization points number in the x and z directions respectively.We omit the detail expressions of other ij n N  n M  A for saving space.Solving (11) by a direct method is impractical when the system is huge.For this reason, various iterative methods, e.g. the Krylov subspace methods, are the interesting alternative.However, Krylov subspace methods are not competitive without a good preconditioned [17].Here we adopt the preconditioned Bi-conjugate gradient stabilized (Bi-CGSTAB) to solve (11).The preconditioned is the so-called "shifted Laplace" preconditioner [18].

Nonoverlapping Domain Decomposition
In no overlapping domain decomposition methods, the computational domain where i is the exterior normal to i .However, the conditions ( 17) may cause ill-posed problem.We propose to use the following equivalent boundary conditions: The sub problem on 1 or 2 together with condition ( 18) or ( 19) is now well posed.We describe the general domain decomposition algorithm in the following.The idea is to adjust iteratively the boundary conditions at the interface between sub domains to obtain the transmission conditions of the type (3) and (17).Initialize , for all then iterating for , where n is the iterative number: Solve the following equations for : ) 2) Solve the following interface equations for ( , ) and for ( , ) x z     : where j  and k  denote the th and kth subdomains respectively.
4) Set 1, n n   return to step (1) to start the next iteration until the solutions convergence.

Conclusions
A new wave simulation technique for the elastic wave equation in the frequency domain is investigated.It is based on the no overlapping domain decomposition method.The computational difference schemes and the corresponding algorithm of domain decomposition are presented.The numerical computations both for a homogeneous model and a three-layered model show the effectiveness of our proposed method.This method can be used in the full-waveform inversion.It can sometimes reduce the computational complexity.
Numerical tests are performed on a homogeneous model and a three-layered model.The source function is set to be a point source.First we consider a homogeneous model with 3000 p v  m/s and m/s.The frequency is 100 Hz.The convergent solutions are shown in Figures 2 and 3.In Figure 2, the strategy of two no overlapping sub domains is adopted.

Figure 2 (
a) is the horizontal displacement and Figure 2(b) is vertical component.In Figure 3, the strategy of three nonoverlapping subdomains is used.

Figure 3 Figure 1 .
Figure 1.Sketch map of domain decomposition.The computational domain is decomposed into no overlapping (a) two subdomains and (b) three subdomains.

Figure 2 .
Figure 2. Horizontal (a) and vertical (b) displacements for a homogeneous model.The solutions are obtained based on the nonoverlapping DDM with two subdomains.The frequency in computations is 100Hz.displacementand Figure 3(b) is the vertical component.The waveform is correct through verification.The solutions on the common interface are almost same and we conclude that the solutions convergence.Next we consider an inhomogeneous model with three-layered media shown in Figure 4. From top to bottom in Figure 4, the velocities of compressive wave are 2000 m/s, 3000 m/s and 4500 m/s respectively, and the velocities of shear wave are 1500 m/s, 2000 m/s and 2500 m/s respectively.The solutions are shown in Fig- ure 5. Figure 5(a) is the horizontal displacement and Figure 5(b) is the vertical displacement.The frequency is still 100 Hz.In computations the strategy of two no overlapping subdomains is adopted.We also present the contours of waveform.The contours corresponding to Figures 5(a) and (b) are shown in Figures 6(a) and (b) respectively.Figure 7 are the solutions with a high frequency 200Hz.Figure 7(a) is the horizontal displace-

Figure 7
are the solutions with a high frequency 200Hz.

Figure 7 (
a) is the horizontal displace-ment and Figure 7(b) is the vertical displacement.Comparing Figure 7 with Figure 6, we see that the waves at high frequency have more oscillation which coincides with the law of physics.

Figure 3 .
Figure 3. Horizontal (a) and vertical (b) displacements for a homogeneous model.The solutions are obtained base on the nonoverlapping DDM with three subdomains.The frequency in computations is 100Hz.

Figure 5 .Figure 6 .
Figure 5. Horizontal (a) and vertical (b) displacements for a three-layered model.The solutions are obtained based on the nonoverlapping DDM with two subdomains.The frequency in computations is 100Hz.

Figure 7 .
Figure 7. Waveform contours of horizontal (a) and vertical displacements for a three-layered model.The solutions are obtained based on the nonoverlapping DDM with two subdomains.The frequency in computations is 200Hz.