Numerical solution of the acoustic wave equation using Raviart–Thomas elements

Journal of Computational and Applied Mathematics 206 (2007) 420 – 431 www.elsevier.com/locate/cam

Numerical solution of the acoustic wave equation using Raviart–Thomas elements Eleanor W. Jenkins∗ Department of Mathematical Sciences, Clemson University, Box 340975, Clemson, SC 29634-0975, USA Received 10 December 2005; received in revised form 27 July 2006

Abstract In this paper we discuss the numerical approximation of the displacement form of the acoustic wave equation using mixed finite elements. The mixed formulation allows for approximation of both displacement and pressure at each time step, without the need for post-processing. Lowest-order and next-to-lowest-order Raviart–Thomas elements are used for the spatial discretization, and centered finite differences are used to advance in time. Use of these Raviart–Thomas elements results in a diagonal mass matrix for resolution of pressure, and a mass matrix for the displacement variable that is sparse with a simple structure. Convergence results for a model problem are provided, as are numerical results for a two-dimensional problem with a point source. © 2006 Elsevier B.V. All rights reserved. MSC: 35L05; 65M60 Keywords: Acoustic wave equation; Mixed finite elements

1. Introduction There has been much recent interest in the numerical simulation of waves using continuous finite element methods [2–5,10,14,18,24], mixed methods [12,13,15,16,27], and discontinuous finite element methods [1,17,22]. Much of this work has focused on simulation of the scalar (pressure) wave equation, and a priori error estimates for mixed formulations of the pressure wave equation were provided in [12,13]. Mixed methods for the displacement form of the wave equation were analyzed in [16], where the a priori estimates obtained required less regularity on the displacement solution than previous estimates. Mixed methods for the vector (displacement) wave equation also use the constitutive relationship that defines the pressure field as the second equation in the formulation. Thus, for certain approximation spaces, one can resolve this constitutive relationship locally. The acoustic wave equation is used to model the effects of wave propagation in heterogeneous media. Using the finite element method for its approximation is attractive because of the ability to handle complex discretizations and because error indicators can be developed to aid in adaptive gridding strategies. The authors in [6,7,27] use mixed methods to approximate acoustic wave behavior in the context of fluid–structure interaction problems. In [6,7], the authors used the Raviart–Thomas lowest-order elements to avoid spurrious modes in the discrete problem that have no physical meaning. ∗ Tel.: +1 864 6566907; fax: +1 864 6565230.

E-mail address: [email protected]. 0377-0427/$ - see front matter © 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.cam.2006.08.003

E.W. Jenkins / Journal of Computational and Applied Mathematics 206 (2007) 420 – 431

421

Wave simulations are also used to interpret seismographic field data and predict damage patterns due to earthquakes [25]. In order to develop an appropriate earth model, geophysicists routinely compare wave simulations against historical data, and assumed models are corrected until the match is deemed acceptable. Because the procedure is implemented within an optimization framework, fast and accurate simulations are required. In this paper, we simulate wave propagations using both lowest-order and next-to-lowest order Raviart–Thomas elements. While many authors have used the lowest-order Raviart–Thomas elements in their simulations, the next-tolowest order elements can be used to obtain higher-order approximations without a substantial increase in the cost of the simulation. We provide numerical validation of previously determined error estimates, and we generate numerical results for a point-source test problem. A limiting form of the more general wave equation is studied here. This model problem is defined as follows. Let  be a bounded domain in Rn , n = 2, 3, with boundary j = D ∪ N . Let n be the outward normal vector to j. The general form of the wave equation is utt − ∇ · ˜ = f

in  × (0, T ),

(1)

∇ ·u=0

on D × (0, T ),

(2)

u·n=0

on N × (0, T ),

(3)

u(·, 0) = u0 ut (·, 0) = u1

in , in ,

(4) (5)

where u is the displacement,  is the density, and ˜ is the stress tensor given by the generalized Hooke’s law ˜ = (∇ · u)I˜ + (∇u + (∇u)T ). Here  > 0 and  are the Lamé coefficients characterizing the material. We let f represent a general source term and let u0 and u1 be the initial conditions on displacements and velocities. The limiting case of (1) with  = 0 is referred to as the acoustic wave equation, which is ˜ = f. utt − ∇ · ((∇ · u)I)

(6)

We assume that  and  are bound above and below by the positive constants 0 , 1 , 0 , and 1 , respectively. This vector equation is equivalent to the scalar wave equation if one makes the substitution p = ∇ · u. The mixed method is established by using this relationship, giving the coupled system utt − ∇p = f,

(7)

−1 p = ∇ · u,

(8)

with the appropriate boundary and initial conditions. The effectiveness of the method analyzed in [16] is demonstrated here by providing simulations using both lowestorder and next-to-lowest-order Raviart–Thomas elements on rectangles [21]. These implementations have several advantages. First, approximations for both vector-valued displacement and pressure are computed at each time step, eliminating the post-processing step typically needed to obtain the displacement estimate from the pressure approximation. Second, each of the mass matrices has a simple structure; in particular, the mass matrix for the pressure equation is diagonal, so that no matrix inversion routine is required. Finally, boundary conditions can be placed on either variable, without any theoretical loss of convergence rates. The rest of this article is organized as follows. An overview of the weak formulation of the model problem and the Raviart–Thomas spaces used in this work are provided in Section 3. Numerical results for model problems are given in Section 4. The first example is used to verify convergence rates. The second example demonstrates the effectiveness of the method by approximating a solution to a point-source problem from the literature. Concluding remarks are provided in Section 5. 2. Notation We use the following inner products and norms in this paper. The L2 inner product over  is defined as  uv d, (u, v) = 

422

E.W. Jenkins / Journal of Computational and Applied Mathematics 206 (2007) 420 – 431

and we denote by  · L2 the L2 norm over , i.e., uL2 () = (u, u)1/2 . The inner product over the boundary j is denoted as  u, v = uv d j

for u, v ∈

H 1/2+ (),

with  > 0. The time–space norm  · L2 (0,T ;L2 ()) is defined as 

uL2 (0,T ;L2 ()) = uL2 (L2 ) =

T 0

1/2 u2L2 ()

.

The time–space norm  · L∞ (L2 ) is similarly defined. In addition to L2 , we also use the standard Sobolev space H(, div) = {v : v ∈ (L2 ())n , ∇ · v ∈ L2 ()}, with associated norm vH(,div) = vL2 + ∇ · vL2 . 3. Weak formulation The weak formulation is derived in the usual manner. In [16], either p = 0 or u · n = 0 is assumed on j, giving the continuous weak formulation (utt , v) + (p, ∇ · v) = (f, v) ∀ v ∈ V, (−1 p, w) − (∇ · u, w) = 0

(9)

∀w ∈ W,

(10)

where V and W represent the standard Sobolev spaces for mixed methods: H(, div) = {v : v ∈ (L2 ())n , ∇ · v ∈ L2 ()}, V = {v ∈ H(, div) : v · n|N = 0}, W = H 1/2+ () for any  > 0. The error estimates in [16] are valid for any of the usual mixed finite element approximating spaces, e.g., Raviart–Thomas–Nedelec spaces [19,21], Brezzi–Douglas–Marini spaces [9], or Brezzi–Douglas–Fortin–Marini spaces [8]. Choosing one of these finite-dimensional spaces as Vh ⊂ V and Wh ⊂ W , where the projection operator for V is defined as h : H(, div) → Vh , and given a time step size t > 0, the fully discrete problem may be stated as follows: find (Un+1 , P n+1 ) in Vh × Wh such that (U0 , v) = ( h u0 , v) ∀ v ∈ Vh ,

(11)

(P , w) = (p , w) ∀ w ∈ Wh ,     2 2  jt U1/2 , v + (P 0 , ∇ · v) = f 0 +  h u1 , v t t

(12)

0

0

∀ v ∈ Vh ,

(13)

(j2t Un , v) + (P n , ∇ · v) = (f n , v) ∀ v ∈ Vh ,

(14)

(−1 P n+1 , w) − (∇ · Un+1 , w) = 0

(15)

∀ w ∈ Wh .

We use the discrete notation jt 1/2 =

1 − 0 , t

with n = (tn ), ti = it.

j2t n =

n+1 − 2 n + n−1 (t)2

,

E.W. Jenkins / Journal of Computational and Applied Mathematics 206 (2007) 420 – 431

423

The Raviart–Thomas spaces on rectangles can be easily described using notation from [23]. For simplicity, let  be the rectangle [0, 1] × [0, 1], and partition  using x : 0 = x0 < x1 < · · · < xn = 1, y : 0 = y0 < y1 < · · · < ym = 1. Define the piecewise-polynomial space Mrq () = {v ∈ C q ([0, 1]) : v is a polynomial of degree r on each subinterval of }, where discontinuous functions are represented by q = −1. The Raviart–Thomas spaces of index r can thus be described using tensor products of these piecewise-polynomial spaces as Whr = Mr−1 (x ) ⊗ Mr−1 (y ), ˜ hr = [Mr+1 (x ) ⊗ Mr−1 (y )] × [Mr−1 (x ) ⊗ Mr+1 (y )], V 0 0 ˜ hr : v1 (0, y) = v1 (1, y) = 0, v2 (x, 0) = v2 (x, 1) = 0} Vhr = {v = (v1 , v2 ) ∈ V ˜ hr : v · n|j = 0}. = {v ∈ V Note that the discrete pressure space Whr is discontinuous, the v1 component of the discrete displacement space is discontinuous at edges parallel to the horizontal axis, and the v2 component is discontinuous at edges parallel to the vertical axis. The bases for M10 (x ) and M10 (y ) are standard piecewise continuous linear functions. Denote these functions by y x vi and vj , respectively. The bases for M0−1 (x ) and M0−1 (y ) are constant functions across the element; denote these y by wix and wj , respectively. The mass matrices for the resolution of the displacement variable are given by y

y

y

y

(Mx )i  j  ,ij = (vix wj , vix wj  ), (My )i  j  ,ij = (wix vj , wix vj  ) and the mass matrix for resolution of the pressure variable is y

y

(Mp )i  j  ,ij = (wix wj , wix wj  ) If the unknowns are ordered by first numbering all of the horizontal unknowns followed by the vertical unknowns, the mass matrices Mx and My are tridiagonal [23,28]. The matrix Mp is diagonal (see [10,11,26]), which means that the pressure unknowns can be resolved by simple division. The basis functions for RT1 are slightly different. Use of the higher-order elements increases the number of unknowns by adding new nodal values. The basis functions M20 (x ) and M20 (y ) are standard piecewise continuous quadratic functions. (Mesh points are now at xi , xi+1 , and xi+1/2 = 21 (xi + xi+1 ) for x , and yj , yj +1 , and yj +1/2 = 21 (yi + yi+1 ) for y .) The basis functions for M1−1 (x ) and M1−1 (y ) are discontinuous piecewise linear functions, and the nodal points for these functions are located at the Gauss points in the element. Sketches of RT0 and RT1 elements are given in Fig. 1 to illuminate this point. (See also [26] for more details.)

Fig. 1. RT0 (left) and RT1 (right) rectangular elements. Nodes for u1 are denoted by ; for u2 , by , and for p by +.

424

E.W. Jenkins / Journal of Computational and Applied Mathematics 206 (2007) 420 – 431 y

Again denoting the basis functions for M1−1 (x ) and M1−1 (y ) by wix and wj , respectively, we note that (wix , wix )=0 y y for i  = i  , and similarly (wj , wj  ), j  = j  . Thus, the mass matrix for pressure is again diagonal (see [11,10,26]). This property also leads to mass matrices Mx and My with a maximum bandwidth of 6 on any row. 4. Numerical results The numerical results provided in this section verify predicted a priori error estimates and show the applicability of the method for a test problem with a point source. 4.1. Convergence rates The error estimates established in [16] were as follows. Theorem 1. If u ∈ L∞ (H(; div)), uttt ∈ L1 (L2 ()) utttt ∈ L∞ (L2 ()), and p ∈ L∞ (L2 ()), then for {Un , P n } 1/2 1/2 defined by (12)–(15) there exists a constant C independent of h and t such that if t < 2h0 /1 C0 , then 1/2 (u − U)l ∞ (L2 ) + 1/2 (p − P )l ∞ (L2 ) C(hr + t 2 )(uL∞ (H r ) + uttt L∞ (L2 ) + pL∞ (L2 ) ),

(16)

where r is associated with the degree of the finite element polynomial. This error estimate can be easily adjusted to consider RTk as the approximation space. Corollary 1. If the approximation space is RTk , then there exists a constant C independent of h and t such that if 1/2 1/2 t < 2h0 /1 C0 , then 1/2 (u − U)l ∞ (L2 ) + 1/2 (p − P )l ∞ (L2 ) C(hk+1 + t 2 )(uL∞ (H k ) + uttt L∞ (L2 ) + pL∞ (L2 ) ).

(17)

Proof. The error estimates found in [16] use approximation properties for the mixed spaces. For the space RTk , we have that [21]  h u − uH(,div) hk+1 (uL2 + ∇ · uL2 ) and Ph p − pL2 hk+1 pL2 , where Ph is the L2 projection of W onto Wh . Thus, the temporal error bounds remain the same, and the spatial error bounds depend on the degree of the lowest-order polynomial associated with the RTk space.  4.1.1. Numerical verification of error estimates The estimates above indicate that the order of convergence of u should be at least that of p. We use the following example to verify the spatial error estimates associated with using RT0 and RT1 for approximating the acoustic wave equation. Example 1. For  = [−1, 1] × [−1, 1], we choose the known solution   2 x −1 , p = ∇ · u = 2x + 2y. u= y2 − 1 Thus u · n = 0 on j.

(18)

E.W. Jenkins / Journal of Computational and Applied Mathematics 206 (2007) 420 – 431

425

Table 1 Approximate errors at t = 1 using u and p and RT0 elements h

u − Uh L2

1 8

p − Ph L2

3.2837e − 2

1 16

8.1408e − 3

2.012

2.0471e − 1

1.005

1 32

2.0247e − 3

2.007

1.0224e − 1

1.002

1 64

5.0010e − 4

2.017

5.1077e − 2

1.001

1 128

1.2027e − 4

2.056

2.5529e − 2

1.001

Rate

Rate

4.1085e − 1

L2 errors vs. mesh size 100 Displacement Pressure

L2 Errors

10-1

10-2

10-3

10-4 10-3

10-2

10-1

100

Mesh Size Fig. 2. L2 errors for u and p using RT0 elements.

The errors listed in Table 1 were obtained at time t = 1, using t = 0.00001; they are shown graphically in Fig. 2. The convergence rates are defined as ln(u − U2h L2 /u − Uh L2 ) ln 2

and

ln(p − P2h L2 /p − Ph L2 ) , ln 2

respectively. These data indicate that the spatial error estimates agree with predicted rates. Notice that the convergence rate for u is slightly higher than what the theory states. The theoretical convergence rates will be dominated by the lowest-order rate, which is that for p. However, since we have that ∇ · u = p, we expect that the convergence rate for u should be one order higher than that for p in order to maintain the theoretical estimate. A similar example is used to verify the error estimate for RT1 . We keep  = [−1, 1] × [−1, 1], and adjust the known solution u as   4 x −1 (19) , p = ∇ · u = 4x 3 + 4y 3 . u= 4 y −1 The errors for RT1 are listed in Table 2, again obtained at t = 1 using t = 0.00001. The errors are plotted in Fig. 3. These values also indicate agreement with predicted convergence rates.

426

E.W. Jenkins / Journal of Computational and Applied Mathematics 206 (2007) 420 – 431

Table 2 Approximate errors at t = 1 using (19) for u and p and RT1 elements h

u − Uh L2

1 8

p − Ph L2

3.4539e − 3

1 16

4.3865e − 3

7.87

2.2801e − 2

3.99

1 32

5.5083e − 5

7.96

5.7053e − 3

4.00

1 64

6.9025e − 6

7.98

1.4267e − 3

4.00

1 128

9.3606e − 7

7.37

3.5661e − 4

4.00

Rate

Rate

9.1001e − 2

L2 errors vs. mesh size 10-1 Displacement Pressure

10-2

L2 Errors

10-3

10-4

10-5

10-6

10-7 10-3

10-2

10-1

100

Mesh Size Fig. 3. L2 errors for u and p using RT1 elements.

4.2. Symmetric point source Example 2. In this example, we mimic a test problem described in [10]. We define =[0, 12]×[0, 12], and approximate the solution of j2 u(x, t) − ∇(∇ · u(x, t)) = g(x)f (t), jt 2 u(x, 0) =

j u(x, 0) = 0, jt

u(x, t) = 0,

x ∈ , t ∈ (0, T ],

x ∈ ,

(20) (21)

x ∈ j, t ∈ (0, T ].

(22)

As in [10], we use f (t) = 2a(2a(t − b)2 − 1) exp(−a(t − b)2 ), where a=

 2 , 1.31

b = 1.35.

t ∈ [0, T ],

(23)

E.W. Jenkins / Journal of Computational and Applied Mathematics 206 (2007) 420 – 431

427

0.04 0.05

0.03 0.02

0

0.01 0

-0.05

-0.01

12

-0.02

10

-0.03 15

8 10 5 0 0

2

4

6

8

10

12

6 4 2 0

2

4

6

8

10

12

Fig. 4. Solution of p for RT0 .

For our simulation, we adjust g(x) slightly, so that ∇ · g(x) is symmetric: 

⎤ −(x − 6) exp −7 (x − 6)2 + (y − 6)2 ⎢ ⎥ g(x) = ⎢ 

⎥ ⎣ ⎦. 2 2 −(y − 6) exp −7 (x − 6) + (y − 6) ⎡

(24)

In [10], the authors solved the scalar wave equation. This is equivalent to the acoustic wave equation under consideration in this work if we take the divergence of both sides of our equation. Therefore, in order for the forcing functions to be similar, the divergence of the forcing function in this work should be as close as possible to the forcing function in [10]. At a minimum, the symmetry of the forcing function should be preserved so that at least visual comparisons with the solutions generated in [10] can be made. We provide graphical solutions to the system described above in the following subsections. Pressure surfaces, displacement surfaces, and values of pressure at a given point in the domain over an extended time period are shown. These solutions were compared visually to the solutions in [10] to assess whether the mixed method is reasonable for the acoustic wave equation. 4.2.1. Results with RT0 We compute the solution for T = 5 using a 200 × 200 mesh with uniform grid spacing in both the x and y directions. The time step is t = 0.001. In Fig. 4, two views of the computed pressure solution p are shown at T = 5. The shape of the pressure wave is consistent with that presented in [10]. The displacement solutions u1 and u2 at T = 5 are given in Fig. 5. The value of p at the point (3, 9) for t ∈ [0, 25] is plotted in Fig. 6 for further comparison with the results in [10]. The behavior of p is similar to that in [10], and the amplitude of the pressure wave differs because of the adjustment to g(x). 4.2.2. Results with RT1 Similar results are obtained using RT1 . In this case, we use a 100 × 100 mesh with a time step t = 0.001. The simulation results at T = 5 for p and u are given in Figs. 7 and 8. The value of p at (3, 9) is plotted for t ∈ [0, 25] (Fig. 9). The amplitude and behavior of the pressure wave remains similar to that obtained in the RT0 simulation, but slight differences begin to emerge as t increases (see Fig. 10). The overall behavior of the wave is also smoother because of the increased order of p.

428

E.W. Jenkins / Journal of Computational and Applied Mathematics 206 (2007) 420 – 431

x 10-3

x 10-3

8 6 4 2 0 -2 -4 -6 -8 15

8 6 4 2 0 -2 -4 -6 -8 15 10 5

6

4

2

0 0

10

8

12

10 5 0 0

6

4

2

8

10

12

Fig. 5. Solution of u = (u1 , u2 ) at T = 5 using RT0 . u1 is on the left, and u2 is on the right.

0.04 0.03

Pressure

0.02 0.01 0 -0.01 -0.02 -0.03 -0.04 0

5

10

15

20

25

Time Fig. 6. Value of p at (3, 9) and t ∈ [0, 25] using RT0 .

0.04

0.05

0.03 0.02

0

0.01 0

-0.05

-0.01

12

-0.02

10

-0.03 15

8 10 5 0 0

2

4

6

8

10

12

6 4 2 0

Fig. 7. Solution of p using RT1 .

2

4

6

8

10

12

E.W. Jenkins / Journal of Computational and Applied Mathematics 206 (2007) 420 – 431

x 10-3

429

x 10-3

8 6 4 2 0 -2 -4 -6 -8 15

8 6 4 2 0 -2 -4 -6 -8 15 10 5

2

0 0

4

8

6

10

12

10 5 0 0

2

Fig. 8. Solution of u = (u1 , u2 ) at T = 5 using RT1 . u1 is on the left, and u2 is on the right.

0.04 0.03

Pressure

0.02 0.01 0 -0.01 -0.02 -0.03 0

5

10

15

20

25

20

25

Time Fig. 9. Value of p at (3, 9) and t ∈ [0, 25] using RT1 .

0.04 0.03

Pressure

0.02 0.01 0 -0.01 -0.02 RT0 RT1

-0.03 -0.04 0

5

10

15 Time

Fig. 10. Comparison between RT0 and RT1 for p at (3, 9) and t ∈ [0, 25].

4

6

8

10

12

430

E.W. Jenkins / Journal of Computational and Applied Mathematics 206 (2007) 420 – 431

Table 3 Conjugate gradient iterations Mesh size

RT space

Total CG Its.

Avg. per time step

50 × 50

RT0 RT1

723,089 795,409

14.5 16

200 × 200

RT0 RT1

489,686 677,281

9.8 13.5

4.2.3. Using conjugate gradient solvers The numerical results presented in the previous sections were obtained using a sparse direct solver for the mass matrix associated with the displacement variable. Since the mass matrix for the displacement variable is symmetric, the conjugate gradient iterative solver can be used. The number of iterations required for the conjugate gradient algorithm to converge to an absolute tolerance of 10−12 are given in Table 3 for a 50 × 50 mesh, using both RT0 and RT1 , and also for a 200 × 200 mesh using both discretizations. We include these numbers to indicate that there is not a substantial amount of increased work in the solution of the linear system, if the cost of a matrix–vector multiplication is assumed to be approximately the same in each case. The time step size for both problems was 0.0001, and there were 50, 000 steps taken in the simulation. The test problem was the point-source test problem. No pre-conditioners were used in either case. 5. Conclusions The error estimates established in [16] provided optimal convergence rates for the use of mixed finite elements methods in solving the acoustic wave equation. The work in this paper numerically verifies those estimates and gives results for a typical point-source problem in a square domain using lowest and next-to-lowest order Raviart–Thomas elements (RT0 and RT1 , respectively) on quadrilaterals. The reduced requirement on the regularity of the displacement variable was a main feature in the original work, as was the ability to resolve both pressure and displacement with no post-processing. Indeed, the pressure variable is determined using the constitutive relationship between displacement and pressure. Additionally, when using RT0 and RT1 spaces for the discretization, the pressure solution requires inversion of a diagonal mass matrix. The diagonal mass matrix for pressure is not unique to the Raviart–Thomas elements, and can be realized by any elements whose nodal points are appropriately defined [10,11]. The results from the point-source test problem indicate that the approach taken in this work is reasonable; indeed, given the simplicity of the mass matrices, it can be expected that these methods will provide increased efficiency in the optimization framework used by geophysical modelers. The model problems also verify the theoretical convergence estimates with regard to both displacement and pressure solutions for the Raviart–Thomas elements. Finally, while a direct solver was used for the simulations presented in this paper, preliminary results using a conjugate gradient iterative solver for the linear system for displacement indicate that only a few additional iterations per time step are required for convergence of the RT1 solution as compared to the RT0 solution. Therefore, if the cost of storing the additional unknowns can be tolerated, and if the cost of a matrix–vector multiplication in both cases is relatively low, it is reasonable to consider using RT1 over RT0 . One may also take advantage of the work done on optimal pre-conditioners for Raviart–Thomas elements, such as those presented in [20]. References [1] M. Ainsworth, P. Monk, W. Muniz, Dispersive and dissipative properties of discontinuous Galerkin finite element methods for the second-order wave equation, J. Sci. Comput. 27 (1–3) (2006) 5–40. [2] G.A. Baker, Error estimates for finite element methods for second order hyperbolic equations, SIAM J. Numer. Anal. 13 (4) (1976) 564–576. [3] W. Bangerth, R. Rannacher, Adaptive finite element techniques for the acoustic wave equation, J. Comput. Acoustics 1 (1) (1999) 1–17. [4] W. Bangerth, R. Rannacher, Finite element approximation of the acoustic wave equation: Error control and mesh adaption, East–West J. Numer. Math. 7 (4) (1999) 263–282.

E.W. Jenkins / Journal of Computational and Applied Mathematics 206 (2007) 420 – 431

431

[5] H. Bao, J. Bielak, O. Ghattas, L.F. Kallivokas, D.R. O’Hallaron, J.R. Shewchuk, J. Xu, Large-scale simulation of elastic wave propagation in heterogeneous media on parallel computers, Comput. Methods Appl. Mech. Engng. 152 (1998) 85–102. [6] A. Bermúdez, R. Durán, M.A. Muschietti, R. Rodríguez, J. Solomin, Finite element vibration analysis of fluid–solid systems without spurious modes, SIAM J. Numer. Anal. 32 (4) (1995) 1280–1295. [7] A. Bermúdez, J.L. Ferrín, A. Prieto, A finite element solution of acoustic propagation in rigid porous media, Internat. J. Numer. Methods Engng. 62 (2005) 1295–1314. [8] F. Brezzi, J. Douglas Jr., M. Fortin, L.D. Marini, Efficient rectangular mixed finite elements in two and three space variables, RAIRO Modèl. Math. Anal. Numér. 21 (1987) 581–604. [9] F. Brezzi, J. Douglas Jr., L.D. Marini, Two families of mixed elements for second order elliptic problems, Numer. Math. 88 (1985) 217–235. [10] G. Cohen, P. Joly, J.E. Roberts, N. Tordjman, Higher order triangular finite elements with mass lumping for the wave equation, SIAM J. Numer. Anal. 38 (6) (2001) 2047–2078. [11] G. Cohen, P. Monk, Gauss point mass lumping schemes for Maxwell’s equations, Numer. Methods Partial Differential Equations 14 (1998) 63–88. [12] L.C. Cowsar, T.F. Dupont, M.F. Wheeler, A priori estimates for mixed finite element methods for the wave equation, Comput. Methods Appl. Mech. Engng. 82 (1990) 205–222. [13] L.C. Cowsar, T.F. Dupont, M.F. Wheeler, A priori estimates for mixed finite element approximations of second-order hyperbolic equations with absorbing boundary conditions, SIAM J. Numer. Anal. 33 (2) (1996) 492–504. [14] T.F. Dupont, L2 -estimates for Galerkin methods for second order hyperbolic equations, SIAM J. Numer. Anal. 10 (1973) 880–889. [15] R. Glowinski, S. Lapin, Solution of a wave equation by a mixed finite element-fictitious domain method, Comput. Methods Appl. Math. 4 (4) (2004) 431–444. [16] E.W. Jenkins, B. Riviére, M.F. Wheeler, A priori error estimates for mixed finite element approximations of the acoustic wave equation, SIAM J. Numer. Anal. 40 (5) (2002) 1698–1715. [17] C. Johnson, Discontinuous Galerkin finite element methods for second order hyperbolic problems, Comput. Methods Appl. Mech. Engng. 107 (1993) 117–129. [18] K.J. Marfurt, Accuracy of finite-difference and finite-element modeling of the scalar and elastic wave equations, Geophysics 49 (5) (1984) 533–549. [19] J.C. Nedelec, Mixed finite elements in R3 , Numer. Math. 35 (1980) 315–341. [20] C.E. Powell, D. Silvester, Optimal preconditioning for Raviart–Thomas mixed formulation of second-order elliptic problems, SIAM J. Matrix Anal. Appl. 25 (3) (2004) 718–738. [21] R.A. Raviart, J.M. Thomas, Mathematical Aspects of the Finite Element Method, Lecture Notes in Mathematics, vol. 106, Springer, Berlin, 1997, pp. 292–315. [22] B. Rivière, M.F. Wheeler, Discontinuous finite element methods for acoustic and elastic wave problems, Contemp. Math. 329 (2003) 271–282. [23] T.F. Russell, M.F. Wheeler, Finite element and finite difference methods for continuous flows in porous media, in: R.E. Ewing (Ed.), Mathematics of Reservoir Simulation, SIAM Publications, Philadelphia, 1983, pp. 35–106. [24] F.J. Sabadell, F.J. Serón, J. Badal, A parallel laboratory for simulation and visualization of seismic wavefields, Geophys. Prospecting 48 (2000) 377–398. [25] M.K. Sen, P.L. Stoffa, Global Optimization Methods in Geophysical Inversion, Elsevier Science, Amsterdam, The Netherlands, 1995. [26] J. Shen, A block finite difference scheme for second-order elliptic problems with discontinuous coefficients, SIAM J. Numer. Anal. 33 (2) (1996) 686–706. [27] X. Wang, K.-J. Bathe, Displacement/pressure based mixed finite element formulations for acoustic fluid–structure interaction problems, Internat. J. Numer. Methods Engng. 40 (1997) 2001–2017. [28] M.F. Wheeler, R. Gonzalez, Mixed finite element methods for petroleum reservoir engineering problems, In: Computational Methods in Engineering and Applied Sciences VI, North-Holland, Amsterdam, 1984, pp. 639–658.