A dispersion minimizing scheme for the 3-D Helmholtz equation based on ray theory
AA DISPERSION MINIMIZING SCHEME FOR THE 3-D HELMHOLTZEQUATION BASED ON RAY THEORY
CHRISTIAAN C. STOLK
Abstract.
We develop a new dispersion minimizing compact finite difference scheme forthe Helmholtz equation in 2 and 3 dimensions. The scheme is based on a newly developedray theory for difference equations. A discrete Helmholtz operator and a discrete operatorto be applied to the source and the wavefields are constructed. Their coefficients arepiecewise polynomial functions of hk , chosen such that phase and amplitude errors areminimal. The phase errors of the scheme are very small, approximately as small asthose of the 2-D quasi-stabilized FEM method and substantially smaller than those ofalternatives in 3-D, assuming the same number of gridpoints per wavelength is used. Innumerical experiments, accurate solutions are obtained in constant and smoothly varyingmedia using meshes with only five to six points per wavelength and wave propagationover hundreds of wavelengths. When used as a coarse level discretization in a multigridmethod the scheme can even be used with downto three points per wavelength. Testson 3-D examples with up to 10 degrees of freedom show that with a recently developedhybrid solver, the use of coarser meshes can lead to corresponding savings in computationtime, resulting in good simulation times compared to the literature. Introduction
We consider the discretization on regular meshes of the Helmholtz equation(1) − ∆ u − k ( x ) u = f with large and variable k . These methods are widely used for simulations on unboundeddomains, for example in exploration geophysics, using domain sizes, in three dimensions,of up to hundreds of wavelengths [23, 5, 22, 25].A key issue for such discretizations are the dispersion (phase) errors, that are closelyrelated to pollution errors [3]. Typically, the propagating wave solutions to the discreteand continuous equations have slightly different wavelengths. These wavelength errors arealso referred to as phase velocity or phase slowness errors, in which case they are differentlynormalized. They lead to phase errors in the solution that grow with the distance fromthe source. A second important consideration is solver cost. The discretized Helmholtzoperator should of course be cheap to apply and/or invert.A class of discretizations, that performs relatively well on these criteria, is given by socalled compact finite difference methods, that use a 3 × × × k and the grid spacing h [3, 13, 16, 21, 27, 6, 29, 26].We will discuss these schemes more in detail below, and compare their phase slownesserrors with those of standard finite differences and Lagrange finite element methods onregular meshes.To design such methods, several strategies have been followed. One approach is tooconstruct schemes of higher order, for example order four order six, see [13, 27, 29] and the Korteweg-de Vries Institute for Mathematics, Science Park 105-107, 1090 XG AmsterdamThe Netherlands
E-mail address : [email protected] . a r X i v : . [ m a t h . NA ] F e b A DISPERSION MINIMIZING SCHEME FOR THE 3-D HELMHOLTZ EQUATION references in [27]. Another approach is to stay with second order schemes but minimize thedispersion errors, because these are the dominant errors for long distance wave propagation[3, 16, 21, 6, 26]. From the point of view of phase slowness errors, the sixth order schemesof [27] and [29] and the quasi-stabilized FEM (QS-FEM) scheme of [3] are the best, see theresults below. The latter has the smallest phase slowness errors by a substantial margin,but is only available in 2-D.Alternative methods include higher order finite elements. An advantage of these meth-ods is the better theory for the behavior of the errors in the limit that the grid spacinggoes to zero, see for example [15, 20].In this paper we introduce a new second order dispersion minimizing scheme in 2 and 3dimensions with phase slowness errors comparable to those of QS-FEM. Accurate ampli-tudes are obtained as well using new amplitude correction operators. A theoretical justifi-cation is given using a newly developed ray theory for Helmholtz-like difference equations.This theory is remarkably similar to the continuous theory, when both are formulated interms of the symbols associated with the operators. With numerical examples we showthe potential for accurate and fast simulation on relatively coarse meshes. In addition weshow applications where the method is used as a coarse level discretization in multigridsolversWe will briefly describe the methodology and the results. It is known that the secondorder, compact finite difference discretizations of the Helmholtz operator form a 3 or 5parameter family, in 2 and 3 dimensions respectively, and that by choosing parameters ina certain way, the phase slowness errors can be reduced compared to standard schemes [16,21]. When coefficients are allowed to depend on hk in a piecewise constant [6] or piecewiselinear fashion [26], they can be further reduced. In this paper we let the parametersdepend in a C fashion on hk through third order Hermite interpolation and obtain afurther reduction of the phase slowness errors.Dispersion minimizing schemes are typically intended for use on quite coarse meshes,and a theoretical understanding that does not involve the limit hk → A ( x ) e iω Φ( x ) . If k is smooth, Φ satisfies a certain eikonal equation and A acertain transport equation, than such solutions approximate the true solutions increasinglywell in the limit ω → ∞ . Here we develop a similar theory for Helmholtz-like differenceequations. We can then choose the discrete scheme such that the phase and amplitudefunctions associated with the discrete operator approximate match those of the continuousoperator well. As can be expected, schemes with small phase slowness errors have accuratephase functions. By introducing amplitude correction operators, accurate amplitude of theray-theoretic solutions are obtained.We are interested in two ways of applying the discretized Helmholtz operators. The firstis simply as a discretization of (1), where the criterion is that the discrete solutions shouldapproximate the true solutions well. Here we are particularly interested in the use of coarsemeshes, say downto five or six points per wavelength, which are for example applied inexploration geophysics [16, 21, 18]. The second application is internally in multigrid basedsolvers. In a multigrid method, the original mesh is coarsened by a factor two one ormore times. On each of the new meshes a discretization of the operator is required. Inthis application the main criterion for a good discretization is that the multigrid methodconverges rapidly. The results concerning the application in multigrid methods are alsoof interest for recently developed two-grid or multigrid methods with inexact coarse levelinverses [5, 25], which are currently some of the fastest solvers in the literature. (Themethod of [25] will actually be tested here.) Below we will write sometimes the fine levelmesh for the original, uncoarsened mesh. DISPERSION MINIMIZING SCHEME FOR THE 3-D HELMHOLTZ EQUATION 3
The small phase slowness errors for IOFD suggest that accurate solutions are possibleeven when quite coarse meshes are used, say downto five or six points per wavelength.We will show that this is indeed the case using numerical examples with constant, andsmoothly varying velocity models (recall that k ( x ) = ωc ( x ) with c the medium velocity).We then consider the application of the IOFD discretization as coarse level discretizationin multigrid based solvers. We will show that in this case IOFD can be used with verycoarse meshes with downto three points per wavelength. With such meshes, solutions aregenerally not accurate enough for direct use, but the approximate solutions can still beused fruitfully in a multigrid method, where they are refined and iteratively improved.This is established using a set of two-dimensional examples, in which a two-grid methodwith IOFD at both levels converges rapidly (see also the results discussed in the nextparagraph). As explained in [26], for the good convergence it is necessary to have very smallphase slowness errors at these very coarse meshes. The IOFD method (in two and threedimensions) and the QS-FEM method (in two dimensions) are the only discretizationsthat have this property to our knowledge, and appear to be uniquely suitable for thisapplication.In 3-D, the fact that a coarser mesh is used does not necessarily imply lower simulationcost. That depends also on the behavior of the solver. To investigate this aspect wepresent tests with a recently developed solver described in [25]. The solver uses a two-gridmethod with an inexact coarse level inverse, given by a double sweep domain decompositionpreconditioner. As described in the previous paragraph, IOFD will also be used as coarselevel discretization. Using the SEG-EAGE Salt Model with up to 10 degrees of freedom asexample, we find that for downto six points per wavelength the cost per degree of freedomchanges little when the frequency is increased. Computation time compare favorably tosome of the results in the literature.The outline of this work is as follows. In section 2 the theory for finite differencediscretizations of the Helmholtz equation with constant k is developed. The symbols andphase slownesses are defined and the discrete Green’s function is studied. In section 3 weconsider the case of variable k and describe ray theory for discrete Helmholtz equations. Insection 4 we compute the phase slowness errors of various existing schemes, as a referencefor the new method. In section 5 we introduce our new interpolated optimized finitedifference method. Section 6 contains some numerical simulations illustrating the accuracyof the solutions when using the IOFD discretization. Section 7 discusses the use of IOFDin multigrid based solvers. Finally, section 8 contains a brief discussion of some furtheraspects. 2. Theory of discrete Helmholtz equations with constant k In this section we study finite difference discretizations of Helmholtz equation(2) Hu = f, H = − ∆ − k in case k is constant. We will assume the grid is given by ( h Z ) d . In this and the nextsection it is convenient to write α, β, ... for multi-indices associated with grid points, suchthat with α = ( α , . . . , α d ) is associated the grid point hα . A difference operator willbe viewed as an operator on functions of x ∈ ( h Z ) d . In this and the next section thedimension d can be any positive integer.For constant k , a finite difference discretization of the Helmholtz operator H is a trans-lation invariant difference operator with coefficients depending on the grid spacing h andon k . By dimensional analysis we may assume that the matrix elements p α,β of such adifference operator P are defined in terms of a finite set of functions f γ by(3) p α,β = 1 h f α − β ( hk ) , A DISPERSION MINIMIZING SCHEME FOR THE 3-D HELMHOLTZ EQUATION where f γ is only nonzero for γ in some finite set Sten( P ).We will first consider the action of such an operator in the Fourier domain and definethe associated symbol and phase slownesses. We next define a “dimensionally reduced”symbol. Then we consider the discrete Green’s function, i.e. solutions to the equation(4) P u = δ, where the δ function on ( h Z ) d is defined by δ ( hα ) = h − d δ α , . . . δ α d , . We obtain thegeneral solution of this equation in the Fourier domain and the asymptotics in the spatialdomain, and determine the same information for the unique outgoing solutions.In the last part of this section we consider a modification of (4) where first a function v is determined that satisfies(5) P v = ˜ Qδ and then u is set equal to(6) u = ˆ Qv.
In this case we assume ˜ Q and ˆ Q are difference operators of order zero. Based on translationinvariance and dimensional reduction, we assume that their matrix elements ˜ q α,β and ˆ q α,β are given by(7) ˜ q α,β = ˜ g α − β ( hk ) , ˆ q α,β = ˆ g α − β ( hk ) , where the ˜ g γ , ˆ g γ are smooth functions that are only nonzero for γ in finite sets Sten( ˜ Q ),Sten( ˆ Q ). A solution u to such a system will be called a modified Green’s function. Ourdiscretization of the Helmholtz equation will be a system of the form (5) and (6), where δ is replaced by the right hand side f .2.1. Symbol and phase slownesses.
To define the symbol, we first define the forwardand inverse Fourier transforms of a function u ( x ), x ∈ ( h Z ) d . They are given by F u ( ξ ) = h d (cid:88) x ∈ ( h Z ) d u ( x ) e − iξ · x (8) F − U ( x ) = (2 π ) − d (cid:90) [ − π/h,π/h ] d U ( ξ ) e iξ · x dξ, (9)where the domain of F u is [ − π/h, π/h ] d . For constant k the finite difference operator P acts like a multiplication in the Fourier domain(10) F ( P u )( ξ ) = P ( ξ ) F u ( ξ ) . where the function P ( ξ ), called the symbol , is given by(11) P ( ξ ) = h − (cid:88) γ f γ ( hk ) e ihγ · ξ . This is similar to the continuous case, where the Helmholtz operator H acts by multipli-cation with H ( ξ ) = ξ − k in the Fourier domain.The Helmholtz equation has propagating plane wave solutions. These are functions u = e ix · ξ that satisfy the homogeneous Helmholtz equation(12) He iξ · x = 0 . They are exactly the plane waves for which ξ is in the zeroset Z H of the symbol of H ( ξ ) = ξ − k (This is of course the set of vectors of length k for H as defined, but theconcept applies more generally.) If P is a translation invariant discretization of − ∆ − k on R d we can similarly look for vectors ξ such that(13) P e ix · ξ = 0 . These are the vectors in the zero set Z P of P ( ξ ). DISPERSION MINIMIZING SCHEME FOR THE 3-D HELMHOLTZ EQUATION 5 If P is a discretization of the Helmholtz operator H then typically the set Z P is closeto, but not identical to Z H . In other words, there are small differences in the wave vectorsof the propagating waves, Z P (cid:54) = Z H . If Z P and Z H can be parameterized by angle, i.e.(14) Z P = { g P ( θ ) θ | θ ∈ S d − } , and similar for Z H (for H ( ξ ) = ξ − k this is of course the case), we define the relativewave number error as a function of θ ∈ S d − by(15) δ ph ( θ ) = g P ( θ ) g H ( θ ) − . Closely related quantities are the phase slowness and phase velocity errors. If k = ωc ( x ) , ξ ∈ Z P , there are associated phase slowness s ph and phase velocity vectors v ph given by(16) s ph = ω − ξ, v ph = ωξ (cid:107) ξ (cid:107) , see e.g. [8]. The quantity δ ph defined in (15) may hence also be called the relative phaseslowness error, or simply phase slowness error.The actual phase error between a numerical and an exact solution is given by (see alsosubsection 2.3)(17) phase error = 2 πδ ph Lλ where L is the distance between source and observation point, λ is the wavelength and δ ph is the phase slowness error associated with the particular angle of propagation. Becauseit is proportional to L/λ the phase error easily may become dominant if is not careful inthe choice of discretization in the high-frequency regime.2.2.
Dimensional reduction and Helmholtz-like symbols.
In case of coefficients,the symbol for arbitrary h can be expressed in terms of that for h = 1(18) P ( ξ ) = 1 h P ( hξ ; hk )where(19) P ( ξ, k ) = (cid:88) γ e iξ · γ f γ ( k ) . By dimensional reduction the symbol ˜ Q ( ξ ), ˆ Q ( ξ ) can be written as(20) ˜ Q ( ξ ) = ˜ Q ( hξ ; hk ( x )) . where(21) ˜ Q ( ξ, k ) = (cid:88) γ e iξ · γ ˜ g γ ( k ) . and similar for ˆ Q .To obtain the results below, we assume that the symbol P is Helmholtz-like as definedin the following
Definition 1.
A symbol P ( ξ ) is said to be Helmholtz-like if the zero set Z P can be param-eterized as in (14), Z P is contained in ] − π/h, π/h [ d , P (0) is negative, ∂P/∂ξ (cid:54) = 0 at allpoints in Z P , and the map (22) N : Z P → S d − : ξ (cid:55)→ ∂P/∂ξ ( ξ ) (cid:107) ∂P/∂ξ ( ξ ) (cid:107) that maps a point in Z P to the unit normal to the surface is a diffeomorphism. It follows that P is Helmholtz-like if P is Helmholtz like. A DISPERSION MINIMIZING SCHEME FOR THE 3-D HELMHOLTZ EQUATION
The discrete Green’s function.
A Green’s function u ( x ) for the discrete equationwill be defined as a solution of (4). If k is constant the equivalent equation for the Fouriertransform U ( ξ ) reads(23) P ( ξ ) U ( ξ ) = 1 . We will first describe the general solution to this equation in the Fourier domain. Thenwe will consider the asymptotics in the spatial domain. Using the results obtained, wecan then derive a unique outgoing Green’s function in the Fourier domain, and state itsasymptotics.Due to the zeros of P , problem (23) has non-unique, distributional solutions. To explaintheir nature, we recall the closely related one dimensional problem to determine all f suchthat(24) xf ( x ) = 1 . We also write this as M x f = 1, where M x is the multiplication operator by the function x . The solutions to (24) are the distributions of the form [12](25) f ( x ) = p . v . x + bδ, where b is a free constant. Here the distribution p . v . x is defined by(26) (cid:104) p . v . x , φ (cid:105) = lim (cid:15) → (cid:90) R \ [ − (cid:15),(cid:15) ] φ ( x ) x dx. In other words δ is in the kernel of M x , while p . v . x is a particular solution to (24).In case of (23) we similarly have a nonzero kernel of M P , with elements BS Z P , where S Z P denotes the singular function of Z P , which is the distribution given by(27) S Z P ( φ ) = (cid:90) Z P φ ( x ) dS ( x ) , and B is any distribution on Z P . For functions f on R d with zero set ˜ Z such that ∇ f ( y ) (cid:54) = 0 for all y ∈ ˜ Z , the principal value u = p . v . f ( y ) can be defined as follows. Let˜ Z (cid:15) = { x | d ( x, ˜ Z ) < (cid:15) } , then(28) (cid:104) p . v . f ( y ) , φ (cid:105) = lim (cid:15) → (cid:90) R d \ ˜ Z (cid:15) φ ( y ) f ( y ) dy. We obtain
Proposition 1.
The solutions to (23) are given by (29) U ( ξ ) = p . v . P ( ξ ) + BS Z P where B can be any distribution on Z P . The freedom in the choice of B is related to the fact that in the spatial domain one canadd any linear combination of plane waves e ix · ξ with ξ ∈ Z P and still have a solution.Let u ( x ) be the inverse Fourier transform of a solution U ( ξ ) for some smooth B (30) u ( x ) = (2 π ) − d (cid:90) [ − π/h,π/h ] d p . v . P ( ξ ) e ix · ξ dξ + (2 π ) − d (cid:90) Z P Be ix · ξ . We will study the asymptotic behavior of this integral for large (cid:107) x (cid:107) using the methodof stationary phase. For p ∈ Z P , we define a certain curvature-like quantity K ( p ) asfollows. After rotating the coordinates, we may assume that ∂P/∂ξ ( p ) is parallel to the d -th coordinate axis and that Z P is locally a graph(31) ξ d = g ( ξ , . . . , ξ d − ) . DISPERSION MINIMIZING SCHEME FOR THE 3-D HELMHOLTZ EQUATION 7
By the assumptions g has a nondegenerate local maximum at ( p , . . . , p d − ). Let − λ , . . . , − λ d − denote the eigenvalues of the second derivative matrix ∂ g∂ ( ξ ,...,ξ d − ) ( p , . . . , p d − ). We de-fine(32) K ( p ) = λ λ . . . λ d − . For d = 3 this is the Gaussian curvature of the surface. In the following proposition N − denotes the inverse function of the map N defined in (22). The result and its proof havesome similarity with results of Lighthill [19]. Proposition 2.
Let u be the inverse Fourier transform of a distribution U as in (29)such that B is a C ∞ function on Z P . Let ξ + = ξ + ( x ) = N − ( x/ (cid:107) x (cid:107) ) and ξ − = ξ − ( x ) = N − ( − x/ (cid:107) x (cid:107) ) . The function u satisfies (33) u ( x ) = (2 π ) − d +12 e − ( d − πi (cid:107) x (cid:107) − d − K ( ξ + ) − / (cid:18) π i (cid:107) ∂P/∂ξ ( ξ + ) (cid:107) + B ( ξ + ) (cid:19) e ix · ξ + + (2 π ) − d +12 e ( d − πi (cid:107) x (cid:107) − d − K ( ξ − ) − / (cid:18) − π i (cid:107) ∂P/∂ξ ( ξ − ) (cid:107) + B ( ξ − ) (cid:19) e ix · ξ − + O ( (cid:107) x (cid:107) − / − d/ ) , (cid:107) x (cid:107) → ∞ . Proof.
We start with the first integral in (30). For x ∈ ( h Z ) d the domain is really a torusand the integrand is C ∞ as a function on the torus. It is convenient to replace the integralon the torus by an integral over a bounded subset of R d . Let ψ be a smooth, positivefunction supported in [ − π/h − η, π/h + η ], that is one on ] − π/h + η, π/h − η [ and satisfies (cid:80) ∞ l = −∞ ψ ( x + 2 πl/h ) = 1 for x ∈ R , and let(34) ψ ( x ) = ψ ( x ) . . . ψ ( x d ) x ∈ R d Then we can write(35) (2 π ) − d (cid:90) [ − π/h,π/h ] d p . v . P ( ξ ) e ix · ξ dξ = (2 π ) − d (cid:90) R d ψ p . v . P per ( ξ ) e ix · ξ dξ for x ∈ ( h Z ) d , where P per is the periodic extension of P , and this formula may also beconsidered for x ∈ R d . We assume η is sufficiently small such that Z P is supported in] − π/h + η, π/h − η [ d .We will write x = τ v , v ∈ S d − and consider the limit τ → ∞ . We assume coordinatesare rotated such that v = (0 , . . . , , χ , denote(36) I χ = (2 π ) − d (cid:90) χ ( ξ ) ψ ( ξ ) p . v . P per ( ξ ) e iτξ d dξ We may assume there are four different types of χ (i) χ = χ + is one on a neighborhood of ξ + (ii) χ = χ − is one on a neighborhood of ξ − (iii) on supp χ ∩ Z P we can write Z P as a graph ξ k = g ( ξ , . . . , ξ k − , ξ k +1 , ξ d )(iv) supp χ ∩ Z P = ∅ We consider these four cases in the limit τ → ∞ using the method of stationary phase[10]. In case (iv) the integral I χ = O ( τ − N ) for any N by the lemma of non-stationaryphase and we don’t need to consider this case further. In case (iii) we can write(37) χψ p . v . P per ( ξ ) = C ( ξ ) p . v . ξ k − g ( ξ , . . . , ξ k − , ξ k +1 , . . . , ξ d ) A DISPERSION MINIMIZING SCHEME FOR THE 3-D HELMHOLTZ EQUATION for some smooth function C ( ξ ) and perform the integral over ξ k . This yields a smoothfunction of ( ξ , . . . , ξ k − , ξ k +1 , ξ d ). By the lemma of non-stationary phase it follows thatagain I χ = O ( τ − N ).In case (i) we can write Z P locally as a graph ξ d = g ( ξ , . . . , ξ d − ). For brevity denote ξ (cid:48) = ( ξ , . . . , ξ d − ) We observe that we can write(38) χψ p . v . P = h ( ξ ) + h ( ξ (cid:48) ) p . v . ξ d − g ( ξ (cid:48) ) − h ( ξ (cid:48) )(1 − ψ ( ξ d − g ( ξ (cid:48) ))) 1 ξ d − g ( ξ (cid:48) ) . where h , h , ψ are smooth, compactly supported functions and ψ = 1 around 0 and h ( ξ (cid:48) + ) = (cid:107) ∂P/∂ξ ( ξ + ) (cid:107) . Then for the first term the lemma of non-stationary phase can beinvoked. Hence this term is O ( τ − N ) for any N . For the third term, the same result canbe obtained using integration by parts. For the second part we recall the standard Fouriertransform F p . v . y = − iπ sgn( η ), it follows that(39) F − p . v . η − a = e iay i y ) . As a consequence, we obtain(40) I χ + = (2 π ) − ( d − i (cid:90) h ( ξ (cid:48) ) e iτg ( ξ (cid:48) ) dξ (cid:48) + O ( τ − N )any N . We can now apply the stationary phase lemma. The function g has its maximumat ξ (cid:48) + and can be expanded as, possibly after a further rotation of coordinates(41) g ( ξ , . . . , ξ d − ) = v · ξ + − d − (cid:88) j =1 λ j ( ξ j − ( ξ + ) j ) + O ( (cid:107) ξ (cid:48) − ξ (cid:48) + (cid:107) ) , see the discussion preceding (32). This yields(42) I χ + = i π ) − ( d − / (cid:107) ∂P/∂ξ ( ξ + ) (cid:107) e − ( d − πi e iτv · ξ + τ − ( d − / K ( ξ + ) − / . The contribution I χ − in case (ii) can be computed similarly, resulting in(43) I χ − = − i π ) − ( d − / (cid:107) ∂P/∂ξ ( ξ − ) (cid:107) e ( d − πi e iτv · ξ − τ − ( d − / K ( ξ − ) − / . For the surface integral(44) (2 π ) − d (cid:90) Z P B ( ξ ) e ix · ξ dS ( ξ ) , we again assume x = τ v and consider the limit τ → ∞ . A partition of unity is applied and,by the method of stationary phase, the only contributions that are not O ( τ − N ) for any N come from neighborhoods of ξ ± . To determine the contribution from a neighborhoodof ξ + , we assume that v is parallel to the d -th coordinate axis so that Z P is locally givenby a graph ξ d = g ( ξ (cid:48) ), ξ (cid:48) = ( ξ , . . . , ξ d − ). The method of stationary phase can be applieddirectly. The only contributions come from neighborhoods of ξ ± , and can be computedsimilarly as for the integral (40).The contribution I χ + , I χ − and the two contributions from the integral (44) togethergive the result. (cid:3) It is straightforward to obtain the outgoing solutions to (4) and (23). In (33) the termwith phase factor e iτv · ξ − must vanish, and we obtain the equation(45) B ( ξ ) = πi (cid:107) ∂P/∂ξ ( ξ ) (cid:107) , for ξ ∈ Z P .We state this as a theorem and include the asymptotic expression for the solution in theresult. DISPERSION MINIMIZING SCHEME FOR THE 3-D HELMHOLTZ EQUATION 9
Theorem 3.
The outgoing solution to (23) is given by (46) U ( ξ ) = p . v . P ( ξ ) + πiS Z P ( ξ ) (cid:107) ∂P/∂ξ ( ξ ) (cid:107) . Its inverse Fourier transform u ( x ) satisfies (47) u ( x ) = (2 π ) − d − e − ( d − πi (cid:107) x (cid:107) − d − i K ( ξ + ) − / (cid:107) ∂P/∂ξ ( ξ + ) (cid:107) e ix · ξ + + O ( (cid:107) x (cid:107) − / − d/ ) , where K and ξ ± are as in proposition 2. The above analysis can be repeated for continuous, Helmholtz like operators with thesame result (see also [19]). For the usual Helmholtz operator H in d = 3 dimensions wehave that Z H given by (cid:107) ξ (cid:107) = k , (cid:107) ∂H/∂ξ (cid:107) | ξ ∈ Z H = 2 k , K + = k and(48) u ( x ) = 14 π (cid:107) x (cid:107) e ik (cid:107) x (cid:107) + O ( (cid:107) x (cid:107) − ) , so that the highest order asymptotic expansion actually equals the well known outgoingGreen’s function.2.4. The modified discrete Green’s function.
Let(49) H ( ξ, k ) = ξ − k . It follows from theorem 3 that if H and P have the same zero sets, i.e. identical phaseslownesses, then the solutions to Hu = δ and P u = δ have asymptotically the same phase.The amplitudes however will differ by a factor (cid:107) ∂H /∂ξ (cid:107)(cid:107) ∂P /∂ξ (cid:107) evaluated at the zero set. In thissubsection we consider therefore the solutions u to the equations (5) and (6), which, aswe will see, obtain different amplitudes.The Fourier transformed solution U ( ξ ) to (5) and (6) is given by the product of thesolution U given in proposition 1 and a factor ˜ Q ( ξ ) ˆ Q ( ξ ). Using this, we can formulate aresult similar to Theorem 3. In this case the adjective outgoing refers to the solution v of(5). The result can be proven by similar arguments as used to prove proposition 2 andtheorem 3. Theorem 4.
The Fourier transform of the outgoing solution to (5) and (6) is given by (50) U ( ξ ) = ˜ Q ( ξ ) ˆ Q ( ξ ) (cid:18) p . v . P ( ξ ) + πiS Z P ( ξ ) (cid:107) ∂P/∂ξ ( ξ ) (cid:107) . (cid:19) . Its inverse Fourier transform u ( x ) satisfies (51) u ( x ) = (2 π ) − d − e − ( d − πi (cid:107) x (cid:107) − d − i ˜ Q ( ξ + ) ˆ Q ( ξ + ) K ( ξ + ) − / (cid:107) ∂P/∂ξ ( ξ + ) (cid:107) e ix · ξ + + O ( (cid:107) x (cid:107) − / − d/ ) , where K and ξ ± are as in proposition 2. Summarizing our findings so far, the discrete solutions u to (5) and (6) are asymptoti-cally equal to the solutions of the continuous Helmholtz equation Hu = δ if the followingtwo conditions are satisfied(i) P ( ξ, k ) and H ( ξ, k ) have the same zero sets(ii) ˆ Q and ˜ Q satisfy(52) ˜ Q ( ξ, k ) ˆ Q ( ξ, k ) = (cid:107) ∂P /∂ξ ( ξ, k ) (cid:107)(cid:107) ∂H /∂ξ ( ξ, k ) (cid:107) for all ( ξ, k ) such that P ( ξ, k ) = H ( ξ, k ) = 0 Theory of discrete Helmholtz equations with variable k In this section we define a class of discrete approximations to the Helmholtz operatorwith variable k , together with the associated symbols. This is the topic of subsection 3.1.We then study ray-theoretic solutions to the equation (4) and to the set of equations (5),(6), where P and Q are now variable coefficient operators.We assume that k = ωc ( x ) , where c is smooth and we consider the limit ω → ∞ . In thediscrete case we assume that hω = constant. Ray-theoretic solutions are then based onthe ansatz(53) u ( x, ω ) = A ( x, ω ) e iω Φ( x ) , for some smoothly varying A and Φ. For the continuous Helmholtz equation, such solutionsare well known and are constructed in two steps. First the ansatz (53) is inserted inthe PDE, and an expansion in ω is performed. Requiring that the highest order termsvanish leads to the eikonal equation for Φ and the transport equation for A . Secondly,initial/boundary conditions for these equations are obtained from the asymptotic behaviorof the constant coefficient solutions. In this way, the solution modulo an error of lowerorder in ω is obtained.For our class of difference equations we follow the same program. The constant co-efficient solutions were already analyzed in subsection 2.3. In subsection 3.2 we find anonlinear first order PDE for Φ and a transport equation for A . Remarkably, we obtainthe same equations in the continuous and discrete case when formulated in terms of thesymbols (which are defined for both continuous and discrete problems). See [19] and [10]for the continuous case and methods used in that case as well as here.In the last part of this section we consider the ray-theoretic solutions to (5) and (6).The conditions (i) and (ii) from subsection 2.4 for P and Q to obtain accurate solutions,need to be modified and extended to have the same ray-theoretic phase and amplitude inthe continuous and discrete case. The operator P should be discretized using a symmetricdiscretization (with = 1 /
2, see below) and we should have ˜ Q = ˆ Q . This is the topic ofsubsection 3.3.3.1. Symbols and operators for variable k . In case k depends on x , finite differencediscretizations of the Helmholtz operator may depend in different ways on the function k . For example, the coefficients p α,β may depend on k and its derivatives at x = hα , butthey may also depend on k at different points, for example on k ( hα ) and k ( hβ ). We willconsider a class of difference operators P , where the matrix elements p α,β depend only onthe value of k at (1 − t ) hα + thβ , where t ∈ { , / , } is a fixed constant. In other wordswe consider operators P with matrix elements of the form(54) p α,β = 1 h f α − β ( hk ((1 − t ) αh + tβh )) . Note that the operator is symmetric if t = 1 / f γ = f − γ . This will turn out to be anappropriate choice for a discrete Helmholtz operator. We will assume that k ( x ) is definedfor all x , not only those in the grid. Similar we assume that for ˜ Q we have(55) ˜ q α,β = ˜ g α − β ( hk ((1 − t ) αh + tβh )) . and similar for ˆ Q .For such operators it is not obvious how to define the symbol. To find an appropriatedefinition, we first consider how to define an operator from a symbol H ( x, ξ ) in the con-tinuous case. This is the subject of pseudodifferential operator theory, and can be donewith the formula [1, 14](56) Op t ( H ( x, ξ )) u = (2 π ) − d (cid:90) (cid:90) H ( x + t ( y − x ) , ξ ) e i ( x − y ) · ξ u ( y ) dξ dy. DISPERSION MINIMIZING SCHEME FOR THE 3-D HELMHOLTZ EQUATION 11
A map from a function H ( x, ξ ) to an operator such as Op t ( H ( x, ξ )) is called a quanti-zation. For t = 0, the previous formula is the standard left-quantization, t = 1 is theright quantization and t = 1 / H ( x, ξ ) = ξ − k ( x ) , thenOp t ( H ( x, ξ )) = − ∆ − k ( x ) , independently of which of these quantizations is used.To obtain a symbol associated with the operator P defined in (54) we rewrite theexpression for P u ( x ) as follows(57) P u ( x ) = h − (cid:88) γ f γ ( hk ( x + tγh )) u ( x + hγ )= h − d (cid:88) γ (cid:88) y ∈ ( h Z ) d f γ ( hk ( x + t ( y − x ))) δ ( x + hγ − y ) u ( y ) . Using the Fourier domain representation δ ( x ) = (2 π ) − d (cid:82) [ − π/h,π/h ] d e ix · ξ this can be rewrit-ten as(58) P u ( x ) = h − d (2 π ) − d (cid:88) y ∈ ( h Z ) d (cid:90) [ − π/h,π/h ] d (cid:88) γ f γ ( hk ( x + t ( y − x ))) e i ( x + hγ − y ) · ξ dξ. This can be written in similar form as (56), namely as(59) Op t ( P ( x, ξ )) u ( x ) def = (2 π ) − d (cid:88) y ∈ ( h Z ) d (cid:90) [ − π/h,π/h ] d P ( x + t ( y − x ) , ξ ) e i ( x − y ) · ξ u ( y ) dξ, where(60) P ( x, ξ ) = h − (cid:88) γ f γ ( hk ( x )) e ihγ · ξ . Thus, associated with P defined in (54) is associated the symbol P ( x, ξ ) given in (60).The parameter t corresponds to the type of quantization, left, right or Weyl quantization.With these definitions, the symbol (60) for variable coefficients may also be expressedentirely in terms of P (61) P ( x, ξ ) = 1 h P ( hξ ; hk ( x )) . A symbol P ( x, ξ ) is called Helmholtz like if it satisfies the definition for each fixed x .3.2. Ray-theoretic equations for amplitude and phase.
In this section we considerthe high-frequency limit ω → ∞ . We assume that ωh = constant, and recall that k ( x ) = ωc ( x ) , where c ( x ) is C ∞ . The operator P and the symbol P ( x, ξ ) become ω -dependent. By˜ P ( x, ξ ) we denote the symbol for ω = 1.(62) ˜ P ( x, ξ ) = 1( ωh ) P ( hωξ, hωc ( x ) ) . For other values of ω we find that(63) P ( x, ξ ; ω ) = ω ˜ P ( x, ξω )We consider the action of P on functions of the form(64) u ( x ) = e iω Φ( x ) A ( x )where Φ and A are C ∞ functions. From the symbol P ( x, ξ ; ω ) and the phase function Φone can derive naturally a vector field, which we call L P, Φ ,ω (cf. [10, section 4.3])(65) (cid:0) L P, Φ ,ω (cid:1) j = ∂P∂ξ j ( x, ω ∇ Φ; ω ) . This vector field is determined by L ˜ P , Φ = ∂ ˜ P∂ξ ( x, ∇ Φ) as follows(66) L P, Φ ,ω ( x ) = ωL ˜ P , Φ ( x ) . Proposition 5.
We have (67) e − iω Φ( x ) P ( e iω Φ( x ) A ( x )) = ω ˜ P ( x, ∇ Φ( x )) A ( x )+ ω i (cid:18) (cid:88) j ( L ˜ P , Φ ) j ∂A∂x j + 12 (div L ˜ P , Φ ) A + ( t − / (cid:88) j ∂ ˜ P∂x j ∂ξ j A (cid:19) + O (1) , ω → ∞ . Proof.
The proof uses a Taylor expansion of the phase function to second order(68) Φ( x + y ) = Φ( x ) + ∇ Φ( x ) · y + 12 (cid:88) j,k ∂ Φ ∂x j x k y j y k + O ( (cid:107) y (cid:107) ) , a Taylor expansion of the amplitude to first order(69) A ( x + y ) = A ( x ) + ∇ A ( x ) · y + O ( (cid:107) y (cid:107) ))and a Taylor expansion of the matrix coefficients to first order(70) p α,β = 1 h f α − β ( hk ( hα )) + th f (cid:48) α − β ( hk ( hα )) h ∇ k ( hα ) · h ( β − α ) . The exponent e iω Φ( x + y ) is then written as a product of three factors(71) e iω Φ( x + y ) = e iω Φ( x ) e iω ∇ Φ( x ) · y iω (cid:88) j,k ∂ Φ ∂x j x k y j y k + O ( ω (cid:107) y (cid:107) ) . These expansions are inserted in the sum(72) (
P u )( x ) = (cid:88) γ h f γ ( hk ( x + thγ )) u ( x + hγ ) . The factor e iω Φ( x ) can be put in front of the expression outside the summation(73) ( P u )( x ) = e iω Φ( x ) h (cid:20) (cid:88) γ e iω ∇ Φ( x ) · hγ f γ ( hk ( x )) A ( x )+ (cid:18) (cid:88) γ e iω ∇ Φ( x ) · hγ f γ ( hk ( x )) ∇ A ( x ) · ( hγ )+ 12 iωA (cid:88) γ e iω ∇ Φ( x ) · hγ (cid:88) j,k ∂ Φ ∂x j ∂x k ( hγ ) j ( hγ ) k f γ ( hk ( x ))+ A (cid:88) γ te iω ∇ Φ( x ) · hγ f (cid:48) γ h ∇ k · ( hγ ) (cid:19) + O ( h ) (cid:21) DISPERSION MINIMIZING SCHEME FOR THE 3-D HELMHOLTZ EQUATION 13
We next use the expression for the symbol P ( x, ξ ) = h (cid:80) γ e ihγ · ξ f γ ( hk ( x )), the followingexpressions for the derivatives of P ( x, ξ ) ∂P∂ξ j = ih (cid:88) γ hγ j e ihγ · ξ f γ ( hk ( x ))(74) ∂ P∂x j ∂ξ k = ih (cid:88) γ ( x γ ) k e ix γ · ξ f (cid:48) γ ( k ( x )) h ∂k∂x j (75) ∂ P∂ξ j ∂ξ k = − h (cid:88) γ ( hγ ) j ( hγ ) k e ihγ · ξ f γ ( hk ( x )) . (76)and an expression for the derivatives of L P, Φ ,ω (77) ∂ ( L P, Φ ,ω ) j ∂x k = ∂ P∂x k ∂ξ j ( x, ω ∇ Φ) + ω (cid:88) l ∂ P∂ξ j ∂ξ l ∂ Φ ∂x l ∂x k . This yields(78) (
P u )( x ) = e iω Φ( x ) (cid:20) P ( x, ω ∇ Φ( x ); ω ) A ( x )+ 1 i (cid:18) n (cid:88) j =1 ( L P, Φ ,ω ) j ∂A∂x j + 12 div L P, Φ ,ω A + ( t −
12 ) (cid:88) j ∂ P∂x j ξ j A (cid:19) + O (1) (cid:21) , ω → ∞ . Using equations (63) and (66) the result follows. (cid:3)
The result is similar to the result in Proposition 4.3.2 of [10]. To find A and Φ suchthat(79) P ( e iω Φ( x ) A ( x )) ≈ P ( x, ∇ Φ) = 0 , which is a nonlinear first order equation like an eikonal equation, and the amplitude mustsatisfy a transport equation(81) (cid:88) j ( L ˜ P , Φ ) j ∂A∂x j + 12 (div L ˜ P , Φ ) A + ( t − / (cid:88) j ∂ ˜ P∂x j ∂ξ j A = 0For t = 1 /
2, this equation conserves | A | .Ray theoretic solutions to equation (4) can now be constructed just as in the continuouscase. By a rescaling, theorem 3 can be used to obtain the asymptotics of a solution for x (cid:54) = 0 and ω → ∞ . The amplitude and phase from formula (47) can hence be used asinitial/boundary values for the eikonal equation for Φ and the transport equation for A ,and these Φ and A can be determined from these equations, where we note that the eikonalequation may not have globally defined solutions, just as in the continuous case.We briefly recall the continuous equivalent of Proposition 5. The following is basically areformulation of proposition 4.3.2 of [10] and can be proven using the method of stationaryphase found in the same text. Proposition 6.
Let H be a continuous Helmholtz like symbol. For the action of Op( H ) on e iω Φ A ( x ) we have the asymptotic development (82) e − iω Φ( x ) Op t ( H )( e iω Φ( x ) A ( x )) = ω ˜ H ( x, ∇ Φ( x )) A ( x )+ ω i (cid:18) (cid:88) j ( L ˜ H, Φ ) j ∂A∂x j + 12 (div L ˜ H, Φ ) A + ( t − / (cid:88) j ∂ ˜ H∂x j ∂ξ j A (cid:19) + O (1) (cid:21) , ω → ∞ . Amplitude correction.
In this section we consider ray-theoretic approximations v = e iω Ψ B and u = e iω Φ A for the solutions v and u to (5) and (6). Assume we havereference ray-theoretic solutions u ref = e iω Φ ref ( x ) A ref ( x ) associated with Helmholtz equa-tion Hu ref = δ where H ( ω ) = − ∆ − ω c ( x ) , i.e. Φ ref satisfies the eikonal equation and A ref the transport equation with appropriate initial conditions. In the following we let H ( ξ, k ) = ξ − k , ˜ H ( x, ξ ) = ξ − c ( x ) .Assume that(i) P ( ξ, k ) and H ( ξ, k ) have the same zero sets(ii) ˜ Q and ˆ Q are identical and Q ( ξ, k ) = ˜ Q ( ξ, k ) = ˆ Q ( ξ, k ) satisfies(83) Q ( ξ, k ) = (cid:107) ∂P /∂ξ ( ξ, k ) (cid:107)(cid:107) ∂H /∂ξ ( ξ, k ) (cid:107) for all ( ξ, k ) such that P ( ξ, k ) = H ( ξ, k ) = 0;(iii) P and H are derived from their respective symbols using t = 1 / We argue that in this case, to highest order u has the same ray-theoretic approximationas u ref . We omit a formal proof, because the arguments are similar as those used above.The construction of the phase and amplitude functions Ψ and B proceeds almost inthe same way as for solutions to (4). Eikonal and transport equations are as follows fromProposition 5. The constant coefficient solutions differ by a factor Q ( ξ + ) and have thesame phase, resulting in different initial/boundary conditions, such that on a small sphereΓ around 0, where we impose the initial/boundary conditions for the eikonal and transportequations, we have(84) Ψ( x ) = Φ ref ( x ) B ( x ) = Q ω =1 ( x, ∇ Φ ref ) A ref ( x )As a result, Ψ( x ) = Φ ref ( x ) everywhere. While we have different transport equations theoperators L ˜ H, Φ ref and L ˜ P , Φ ref are scaled versions of each other(85) L ˜ P , Φ ref = Q ω =1 ( x, ∇ Φ ref ) L ˜ H, Φ ref , It follows from this fact and the transport equation for A ref , that(86) (cid:88) j L ˜ P , Φ ref ∂Q − ω =1 A ref ∂x j + (div L ˜ P , Φ ref ) Q − ω =1 A ref = 0 . This and (84) shows that(87) B ( x ) = Q − ( x, ∇ Φ ref ) A ref ( x )everywhere.The function u is given by applying Q to v . The action of Q on e iω Ψ B is to highestorder equal to a multiplication by Q ( x, ω ∇ Ψ; ω ), so that(88) Φ = Ψ = Φ ref and A ( x ) = ˜ Q ( x, ∇ Φ ref ) B ( x ) = A ref ( x ) , DISPERSION MINIMIZING SCHEME FOR THE 3-D HELMHOLTZ EQUATION 15 concluding the argument.4.
Phase slowness errors for existing discretizations
In this section we will describe three types of discretizations of the Helmholtz equa-tion (1), namely standard finite differences, compact finite differences and Lagrange finiteelements on regular meshes. We then compute phase slowness errors to compare the per-formance of the different methods in this respect, and to obtain reference values for ournew method constructed below.We modify the notation compared to the previous two section. In this section the degreesof freedom for all three types of methods are denoted by u j,k,l (in three dimensions) andassociated with a regular mesh with grid spacing h . For finite element methods of order N the cells are of size N h , and cell boundaries are located at j, k, l ≡ N . Occasionallywe will use d = 2 or 3 to denote the dimension of space.4.1. Standard finite differences.
In a standard finite difference discretization of theoperator − ∆ − k each of the one-dimensional second derivatives in the Laplacian ∆ = ∂ ∂x + ∂ ∂x + ∂ ∂x is approximated by a central difference approximation of the given order.These are given by(89) D ( N )2 u l = h − N/ (cid:88) m = − N/ c ( N ) m u l + m where the c ( N ) m are as in the following table for N = 2 , , , m = − N = 2 1 -2 14 −
112 43 −
52 43 − −
320 32 − −
320 190 − −
15 85 − −
15 8315 − The discrete approximation to the term − k ( x ) u in (1) is simply given by − k l,m,n u l,m,n .The two-dimensional case can be done similarly.4.2. Compact finite difference discretizations.
For constant k , compact finite differ-ence discretizations take the form(90) ( Au ) l,m,n = (cid:88) ( p,q,r ) ∈{− , , } a p,q,r u l + p,m + q,n + r . Because of symmetry, there are four different coefficients A j , j = 0 , , , a p,q,r = A | p | + | q | + | r | , ( p, q, r ) ∈ {− , , } In 2-D we have(92) ( Au ) l,m = (cid:88) ( p,q ) ∈{− , , } a p,q u l + p,m + q . and there are three different coefficients A j , j = 0 , , a p,q = A | p | + | q | , ( p, q ) ∈ {− , , } . The choice of coefficients is done in different ways in [3, 13, 16, 21, 27, 6, 29, 26]. TheQS-FEM method [3] is a two-dimensional method, for which the coefficients are given modulo an overall normalization by(94) A QS − FEM0 = 4 A QS − FEM1 = 2 c ( α ) s ( α ) − c ( α ) s ( α ) c ( α ) s ( α )( c ( α ) + s ( α )) − c ( α ) s ( α )( c ( α ) + s ( α )) A QS − FEM2 = c ( α ) + s ( α ) − c ( α ) − s ( α ) c ( α ) s ( α )( c ( α ) + s ( α )) − c ( α ) s ( α )( c ( α ) + s ( α ))where α = kh and the auxiliary functions c , s , c , s are defined by(95) c ( α ) = cos (cid:16) α cos π (cid:17) s ( α ) = cos (cid:16) α sin π (cid:17) c ( α ) = cos (cid:18) α cos 3 π (cid:19) s ( α ) = cos (cid:18) α sin 3 π (cid:19) . We will not discuss the fourth order method of [13] because it contains still a free parameterand one of the authors has later published a sixth order method in [29]. In the lattermethod variations of k are taken into account. In case of constant k , the coefficients aregiven in three dimensions by(96) A CHO60 = + 6415 − k h
15 + k h A CHO61 = −
715 + k h A CHO62 = − − k h A CHO63 = − A CHO60 = 103 − k h
45 + k h A CHO61 = − − k h A CHO62 = − − k h , again modulo an overall constant. For the method of Sutmann [27] we have (this methodis only for 3-D)(98) A SUT0 = 6415 (cid:0) − k h + k h − k h (cid:1) A SUT1 = − (cid:0) − k h (cid:1) A SUT2 = − (cid:0) k h (cid:1) A SUT3 = − − ∆ − k are split in a contribution from − ∆ anda contribution from − k . We define a three dimensional, symmetric discretization of theidentity M depending on three parameters α , α , α by(99) ( M u ) l,m,n = (cid:88) ( p,q,r ) ∈{− , , } m p,q,r u l + p,m + q,n + r m p,q,r = M | p | + | q | + | r | , ( p, q, r ) ∈ {− , , } where now M = α , M = α , M = α , M = − α − α − α . For constant k the discreteform of the term − k u is given by − k ( M u ) l,m,n . Before defining the negative Laplacianwe define a two-dimensional weighting operator, discretizing the identity, depending ontwo additional parameters α , α (100) ( N [1 , u ) l,m,n = (cid:88) p,q ∈{− , , } n p,q u l + p,m + q,n n p,q = N | p | + | q | , ( p, q ) ∈ {− , , } with N = α , N = α , N = − α − α . Each of the second deratives in the Helmholtzoperator will be discretized using the tensor product of a two-dimensional weighting oper-ator that discretizes the identity and the standard second order discrete second derivative D [ j ]2 , where j = 1 , − D [1]2 ⊗ N [2 , − D [2]2 ⊗ N [1 , − D [3]2 ⊗ N [1 , − k M In the two-dimensional case there are in total three parameters α , α , α , with M = α , M = α , M = − α − α , and N = α , N = − α . The coefficients A j for the 3-D caseare given in terms of the α j (modulo an overall constant h − ) by(102) A = 6 α − ( kh ) α A = − α + (1 − α − α ) − ( kh ) α A = − α + α − ( kh ) α A = − (1 − α − α ) − ( kh ) (1 − α − α − α ) . For the 2-D case we have(103) A = 4 α − ( kh ) α A = 1 − α − ( kh ) α A = − α − ( kh ) (1 − α − α ) . An advantage of this formulation using tensor products is that in case of PML layersaligned with the coordinate axes, the second order operator D [ j ]2 can simply be replacedby its PML-modified version .In this way we have derived a family of second order accurate discretizations. In [16]and [21] the same family of discretizations is considered in two resp. three dimensions(but differently parameterized), and coefficients are chosen such that the maximum phaseslowness error is minimized, where the maximum is taken over all angles and a range of kh corresponding to at least four points per wavelength. This leads to the choices(104) α OPT1 = 0 . α OPT2 = 0 . α OPT3 = 0 . α OPT4 = 0 . α OPT5 = 0 . α JSS1 = 0 . α JSS2 = 0 . α JSS3 = 0 . α j are allowed to vary. In [6] a set of 7 parameters (in three dimensions)is chosen piecewise constant. We will not describe this method in detail but refer to thepaper for resulting phase errors. In [26] the above described set of 5 parameters are chosenas piecewise linear functions. However, in this work, the aim is different, because the phaseslowness differences with a fine scale operator are minimized, not with the exact operator,so that the values of the phase slowness errors cannot be compared.4.3. Lagrange finite elements on regular meshes.
For the description of Lagrangefinite elements on regular meshes, which will also be used in some of the numerical exam-ples, we start with the one-dimensional case. In this case the finite element cells are theintervals (( j − N h, jN h ), j = 1 , , . . . , each containing N − , N ), and shape functions on this reference cell In a PML layer, say a layer associated with x = constant, the derivative ∂∂x is replaced bya α PML , ( x ) ∂∂x where α PML , ( x ) = iσ ( x ) /ω and the function σ indicates the local amount ofdamping [17, 4, 7]. We choose σ quadratically increasing. The discrete second derivative in thefirst coordinate in this PML layer becomes h − α PML , ( x l,m,n )( α PML , ( x i +1 / ,j,k )( u l +1 ,m,n − u l,m,n ) − α PML , ( x i − / ,j,k )( u l,m,n − u l − ,m,n )) By rescaling the equations with a factor α − , the symmetry ofthe system is restored. are given by standard Lagrange polynomials, which we will denote by L ( N ) j ( x ), and whichare one at x = j and zero at x = 0 , , . . . , j − , j + 1 , . . . N − , N and defined to be 0outside [0 , N ]. Letting k ∈ Z and l ∈ { , , . . . , N − } , the one-dimensional trial and testfunctions are given by(106) ψ ( N ) kN + l ( x ) = (cid:40) L ( N )0 ( xh − ( k − N ) + L ( N ) N ( xh − kN ) if l = 0 L ( N ) l ( xh − kN h ) otherwiseThe two and three-dimensional trial and testfunctions are given by tensor products of the ψ ( N ) j . The finite element discretizion is of course derived from the weak form(107) Φ( u, v ) def = (cid:90) Ω d (cid:88) j =1 ∂u∂x j ∂v∂x j dx − (cid:90) Ω k uv dx = (cid:90) Ω f v dx. The elements of the matrix in the discretization are given by(108) a l,m,n ; p,q,r = Φ( ψ p,q,r , ψ l,m,n )The contribution from the term (cid:82) Ω (cid:80) dj =1 ∂u∂x j ∂v∂x j dx can be called the stiffness matrix andthe contribution from (cid:82) Ω k uv dx can be called the mass matrix. If k is constant (or cellwiseconstant) , the stiffness and mass matrices can be computed exactly. If k is variable, thenonly the stiffness matrix can be computed exactly, and for the mass matrix some sort ofquadrature must be used. For constant k these computations are standard and easily doneusing a computer algebra system, and we will not write down the resulting coefficients.For constant k , the finite difference methods are obviously translationally symmetric,i.e. if we denote by a l,m,n ; p,q,r the matrix elements we have(109) a l,m,n ; p,q,r = a l + A,m + B,n + C ; p + A,q + B,r + C For the finite elements there is a symmetry under a subset of translations given by the
A, B, C that are multiples of N .4.4. Phase slowness errors.
For finite difference methods, finding the phase velocitiesor slownesses comes down to determining the zeros of the symbol P ( ξ ) associated witha difference operator P , see subsection 2.1. The symbol is not difficult to obtain, forexample, for compact finite difference discretizations, the symbol is(110) P ( ξ ) = h − (cid:2) A + 2 A (cos( hξ ) + cos( hξ ) + cos( hξ )) + 2 A (cos( h ( ξ + ξ )) + cos( h ( ξ − ξ ))+ cos( h ( ξ + ξ )) + cos( h ( ξ − ξ )) + cos( h ( ξ + ξ )) + cos( h ( ξ − ξ )))+ 2 A (cos( h ( ξ + ξ + ξ )) + cos( h ( ξ − ξ + ξ )) + cos( h ( ξ + ξ − ξ )) + cos( h ( ξ − ξ − ξ ))) (cid:3) To compute the zeros numerically, the standard numerical solver fsolve from Matlab wasused, as well as the more accurate version vpasolve.For finite element methods the elements of the kernel of the operator are no longersimple plane waves, but Bloch waves. In the appendix is described how we compute thephase slowness errors in this case.Phase slowness errors are directionally dependent, i.e. they depend on θ ∈ S d − asexplained in section 2.1. They also depend on kh or equivalently on the number of pointsper wavelength G = πkh . We have computed the maximum relative phase slowness errorsover θ ∈ S d − for a number of schemes as a function of 1 /G . In two and three dimensionsthese schemes are the finite element schemes of order 1 , , , , , , DISPERSION MINIMIZING SCHEME FOR THE 3-D HELMHOLTZ EQUATION 19 r e l a t i v e pha s e s l o w ne ss e rr o r -8 -7 -6 -5 -4 -3 -2 -1 phase slowness errors for various 2-D schemes FE1FE2FE3FE4FE6FE8FD2FD4FD6FD8JSSCHO6QS-FEM
Figure 1.
Phase slowness errors for some 2-D schemes as a function ofthe inverse number of points per wavelength 1 /G . pha s e s l o w ne ss e rr o r -8 -7 -6 -5 -4 -3 -2 -1 phase slowness errors for various 3-D schemes FE1FE2FE3FE4FE6FE8FD2FD4FD6FD8OPT4CHO6SUT
Figure 2.
Phase slowness errors for some 3-D schemes as a function ofthe inverse number of points per wavelength 1 /G .have the method of Operto Et Al [21], indicated by OPT4 and the method of Sutmann[27], indicated by SUT. We have not included results on the method of [6], these are givenin Figure 2(c) in that work. The phase slowness errors as a function of 1 /G are plotted inFigure 1 and 2. At the end of section 5 we will briefly discuss these results.5. A dispersion minimizing scheme with amplitude corrections
In this section we will define our new discretization of the Helmholtz equation. In thisscheme, the approximate solution to the Helmholtz equation Hu = f , H = − ∆ − k ( x ) isfound by solving a discrete system(111) P v = Qf and then setting(112) u = Qv, where P and Q are compact finite difference operators defined momentarily. In section 3 we studied ray-theoretic solutions to difference equations of the type (111)and (112) with f = δ , and we observed that the ray-theoretic solution to these equationwould be identical to those of the Helmholtz equation(113) Hu = δ, H = − ∆ − k ( x ) if the requirements (i) to (iii) of subsection 3.3 are satisfied. It follows from the derivationsthat if these properties are not satisfied exactly, but there are small differences betweenthe zero set of P and that of H and between the values of Q ( ξ, k ) and ∂P /∂ξ ( ξ,k ) ∂H /∂ξ ( ξ,k ) thenthere will be small errors in the phase and amplitude of the ray-theoretic solutions. Theoperators P and Q will be chosen such that these differences are minimal. We will firstconstruct P in subsection 5.1. Then Q will be constructed in subsection 5.2. In subsectionwe will discuss the phase errors of the new method.5.1. IOFD discretization of the Helmholtz operator.
In sections 2 and 3 a generalform for P was given in terms of functions f γ of kh , see (3), (54). For the 5 or 3 dimensionaloperator family of subsection 4.2 (in 3 and 2 dimensions respectively), these are given by(114) f γ = α − ( kh ) α for γ ∈ {− , , } , | γ | = 0 − α + α − ( kh ) α for γ ∈ {− , , } , | γ | = 1 − α + (1 − α − α ) − ( kh ) α for γ ∈ {− , , } , | γ | = 2 − (1 − α − α ) − ( kh ) (1 − α − α − α ) for γ ∈ {− , , } , | γ | = 3,in 3-D and by(115) f γ = α − ( kh ) α for γ ∈ {− , , } , | γ | = 01 − α − ( kh ) α for γ ∈ {− , , } , | γ | = 1 − α − ( kh ) (1 − α − α ) for γ ∈ {− , , } , | γ | = 2.in 2-D. Here | γ | = | γ | + . . . + | γ d | and we used equations (102) and (103). We let α j dependon hk π = 1 /G , where G is the number of points per wavelength used in the discretization(116) α j = α j (1 /G ) , /G = kh π . Next we will choose a parameterization for these function and we will describe how, byminimizing the phase slowness errors in a least-squares sense, we obtain suitable choicesof the functions α j , j = 1 , . . . , d − α j where chosen to depend piecewise linearly on 1 /G . Here we let α j dependpiecewise polynomially on 1 /G , using Hermite interpolation. We will specify a numberof control nodes, and at each node the value of α j and its first derivative ∂α j ∂ (1 /G ) areprescribed. We will assume that the coefficients α j vary slowly, so that we can indeeddefine the four coefficient of the stencil using five parameters depending on 1 /G . If n C denotes the number of control nodes, in this way the functions α j are parameterized by2 n C parameters. In total we have (4 d − n C parameters, collectively denoted by P .Next we specify the objective functional. The first contribution to the objective func-tional is the square integrated phase slowness error, integrated over angle and 1 /G .Because of the symmetries, the phase slowness error need not be integrated over all θ ∈ S d − , but can be integrated over a subset Θ d of the sphere. In 2 dimensions, theangle variable can be chosen in Θ = [0 , π/ θ = ( θ , θ ) = ( polar angle, azimuthal angle), the domain is Θ = [0 , π/ × [0 , π/ dα j d (1 /G ) .In summary, we have(117) T ( P ) = (cid:90) /G max (cid:90) Θ d | δ ph ( θ, P ) | dθ d (1 /G ) + λ (cid:90) /G max d − (cid:88) j =1 (cid:12)(cid:12)(cid:12)(cid:12) dα j d (1 /G ) (cid:12)(cid:12)(cid:12)(cid:12) d (1 /G ) . DISPERSION MINIMIZING SCHEME FOR THE 3-D HELMHOLTZ EQUATION 21 /G α ∂α ∂ (1 /G ) α ∂α ∂ (1 /G ) α ∂α ∂ (1 /G ) Table 1.
Coefficients two-dimensional IOFD /G α ∂α ∂ (1 /G ) α ∂α ∂ (1 /G ) α ∂α ∂ (1 /G ) α ∂α ∂ (1 /G ) α ∂α ∂ (1 /G ) Table 2.
Coefficients three-dimensional IOFDThis integral is discretized, using a weighted sum of regularly sampled contributions. By n A we denote the number of angles to discretize the integration over the sphere and by n G the number of choices for the parameter 1 /G . In two dimensions we chose n A = 20, inthree dimensions n A = 200. The 1 /G axis was discretized in steps of 0 .
01 in the integral(117). The value of λ = 10 − was used and control nodes where chosen in the interval[0 , .
4] with distance 0 . /G in[0 , .
2] were determined, then for 1 /G in [0 . , .
3] keep the values for 1 /G < . /G in [0 . , .
4] keeping those for 1 /G < . G ∈ [0 , . /G ∈ [0 , .
1] were most sensitive to details of themethod. As the phase speed errors are very small in this parameter range we have notexplored this further.The results of the optimization for the two- and three-dimensional case are given inTables 1 and 2. The phase slowness errors as a function of 1 /G (maximum over angle) aregiven in Figure 3, together with those of the IOFD and CHO6 methods.5.2. Amplitude correction operators.
Here we will construct the difference operators Q . We will use the second order discretizations of the identity defined in (99). Here wewill denote the coefficients of this family by β j , which will be functions of kh π = 1 /G . Thus r e l a t i v e pha s e s l o w ne ss e rr o r -8 -6 -4 -2 phase slowness errors of IOFD, QS-FEM, SUT and CHO6 QS-FEMCHO6 2D and 3DSUTIOFD-2DIOFD-3D
Figure 3.
Phase slowness errors for the IOFD method compared to QS-FEM and the sixth order method of [27] and [29].the functions g γ associated with Q , see (7) are given by(118) g γ = β for γ ∈ {− , , } , | γ | = 0 β for γ ∈ {− , , } , | γ | = 1 β for γ ∈ {− , , } , | γ | = 2 − β − β − β for γ ∈ {− , , } , | γ | = 3,in 3-D and by(119) g γ = β for γ ∈ {− , , } , | γ | = 0 β for γ ∈ {− , , } , | γ | = 1 − β − β for γ ∈ {− , , } , | γ | = 2,in 2-D, where β j = β j (1 /G ). The β j (1 /G ) will be defined by Hermite interpolation fromcontrol values similarly as we did for the α j (1 /G ). The control values are chosen tominimize a discrete approximation of the integral(120) (cid:90) /G max (cid:90) Θ d (cid:34) Q ( ξ ) − (cid:115) (cid:107) ∂P/∂ξ ( ξ ) (cid:107)(cid:107) ∂H/∂ξ ( ξ ) (cid:107) (cid:35) ξ = ωs ph ( θ ) θ dθ d (1 /G ) . This integral is discretized in the same way as in the previous subsection. This results ina linearly constrained linear least squares problem which is easy to solve in Matlab. Theresulting coefficients are given in Tables 3 and 4 below. The maximum over angle of theerror Q ( ξ ) (cid:113) (cid:107) ∂P/∂ξ ( ξ ) (cid:107)(cid:107) ∂H/∂ξ ( ξ ) (cid:107) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ξ = ωs ph ( θ ) θ varied between around 10 − for 1 /G = 0 .
05 and 10 − for 1/G= 0.4.5.3. Comparison of phase slowness errors.
The following conclusions can be drawnfrom the data in Figures 1, 2 and 3. First the QS-FEM method of [3] (in two dimensions)and the IOFD method developed here (in two and three dimensions) perform remarkablywell considering their small stencils. They provides a substantial improvement, roughly a
DISPERSION MINIMIZING SCHEME FOR THE 3-D HELMHOLTZ EQUATION 23 /G β ∂β ∂ (1 /G ) β ∂β ∂ (1 /G ) Table 3.
Coefficients amplitude correction operator Q in 2-D /G β ∂β ∂ (1 /G ) β ∂β ∂ (1 /G ) β ∂β ∂ (1 /G ) Table 4.
Coefficients amplitude correction operator Q in 3-Dfactor 20, in phase errors compared to the compact sixth order scheme SUT and CHO6 of[27] and [29], which in turn are better than other alternatives. For higher order FD andFE methods, as can be expected, the error becomes small if both the number of pointsper wavelength and the order N become large, however this effect sets in quite late, e.g.at eight points per wavelength and N = 8 the relative phase slowness errors of the finiteelement method are roughly equal to those of QS-FEM and IOFD.Next we discuss how much accuracy might be needed, and in how far the improvementswill make a difference in simulations. In view of (17) it is not unreasonable to require atleast that δ ph (cid:46) . λL . In a regime of wave propagation over several hundreds wavelengths,using a mesh with five points per wavelength, from the methods considered only QS-FEMand IOFD satisfy this. At six points per wavelength the CHO6 method is near thisbound while FE8 (which is much more expensive) also qualifies. So in these situations theimproved phase slowness accuracy obtained by using QS-FEM or IOFD can be expectedto have some impact in terms of lower cost compared to FE8 and in terms of improvedaccuracy compared to CHO6 and other compact finite difference methods. The latter willbe confirmed in the examples in the next section.6. Numerical examples
In this section we present two numerical experiments, first in a constant medium, andthen in a smoothly varying medium. We will present two-dimensional examples with largedomain sizes on the order of hundreds of wavelengths.As mentioned, phase slowness errors typically lead to phase shift errors in the solutions.Considering wave propagation over 500 wavelengths as an example, it follows from (17) andthe surrounding discussion that these phase shifts errors for IOFD should be negligiblysmall for meshes with five or six points per wavelength, and still quite small for fourand three points per wavelength. For other methods these errors should show up much (a) y (b) theta r Figure 4.
Phase shift errors are easily observed by plotting a 45-degreepart of an annulus, see figure (b). Cubic spline interpolation is applied tomap the data of figure (a) to polar coordinates.stronger. In our first example we will verify this numerically, assuming a constant velocitymodel.To simulate a point source at a given grid point, we will simply use a discrete δ -function.An unbounded domain is simulated by adding a damping layer around the domain ofinterest, with a nonzero imaginary contribution to k that quadratically increases fromthe boundary of the domain of interest In a 1-D damped Helmholtz equation − dudx − k u with k constant, k = α + iβ solutions decay as u = e i ( α + iβ ) x . If k varies slowly the damping becomes proportional to e − (cid:82) β ( x ) dx . The quadratic profileis chosen such that e − (cid:82) β ( x ) dx is on the order of 0 .
001 to 0 .
01. Reflected waves pass twice through thedamping layer. Unfortunately reflections occur due to the medium variations. To make these small Im( k )must increase slowly and these layers must be quite thick. In our experiments we used on the order of 5to 10 wavelengths. DISPERSION MINIMIZING SCHEME FOR THE 3-D HELMHOLTZ EQUATION 25 (a) theta r polar plot FD2 at 10 ppw0 0.2 0.4 0.62020.220.420.620.8 (b) theta r polar plot JSS at 6 ppw (c) theta r polar plot CHO6 at 6 ppw0 0.2 0.4 0.6500500.2500.4500.6500.8 (d) theta r polar plot CHO6 at 5 ppw (e) theta r polar plot CHO6 at 4 ppw0 0.2 0.4 0.6500500.2500.4500.6500.8 (f) theta r polar plot IOFD at 6 ppw (g) theta r polar plot IOFD at 5 ppw (h) theta r polar plot IOFD at 4 ppw (i) theta r polar plot IOFD at 3 ppw (j) theta r polar plot IOFD at 2.5 ppw Figure 5.
Plots of numerical solutions over a 45 degree part of an annulusfor several numerical methods. (a) FD2 at 10 ppw; (b) JSS at 6 ppw; (c),(d), (e) sixth order method of [29] at 6, 5 and 4 ppw; (f)-(j) IOFD methodat 6, 5, 4, 3 and 2.5 ppw.methods). At 6, 5 and ppw the maximum phase errors at 500 wavelengths are 0.27, 0.84and π radians respectively and the associated phase shifts are increasingly visible in thepictures. In parts (e) to (i) we plot results for the IOFD method at 6, 5, 4, 3 and 2.5ppw. At 6, 5 and 4 ppw the maximum phase shifts are 0.0065, 0.020 and 0.089 radiansrespectively, i.e. considerably smaller than observed for CHO6. At 3 ppw the phase shiftafter 500 wavelengths is clearly visible, only at 2.5 points per wavelength does it becomelarge and in this case the field is plotted at 100 instead of less than 500 wavelengths fromthe source.In Figure 6 we plot the amplitude errors for IOFD at 3 and 4 ppw. For more than 4ppw they were increasingly small.In our second example k is variable. To avoid that errors due to the discretizationof the velocity model become dominant we use use a smoothly varying velocity model,namely a smoothed Marmousi model. In this example we will compare a solution withIOFD using a minimum of six points per wavelength with a fourth order finite elementsolution using twice as many grid points in each direction. In these examples the righthand side was a point source and the linear systems were again solved with MUMPS. Incase of variable coefficients we assumed that k ( x ) is defined on the cell centers. The values (a) theta r amplitude errors IOFD at 4 ppw -4 -505 (b) theta r amplitude errors IOFD at 3 ppw -3 -505 Figure 6.
Amplitude errors relative to exact solution for IOFD at (a) 4ppw (b) 3 ppw. smoothed Marmousi velocity modelx(m) y ( m ) Figure 7.
Smoothed Marmousi velocity modelof f γ ( kh ( x )) (see (54)) at other points were obtained using linear interpolation from thevalues of f γ ( hk ( x )) at the cell centers.The velocity model is given in Figure 7. It is obtained from the Marmousi model byconvolving along both of the axes with a cos square pulse of width 160 meter. We willgive results for 50 and 100 Hz. A solution for the first case is given in Figure 8. Figure 9contains four plots. The top plots are reference amplitudes for obtaining relative errorsand the bottom two plots are relative errors with respect to the reference values. In bothcases we give results for 50 and for 100 Hz. The reference value is a local average of theabsolute value of the solution over a square of about 2 by 2 wavelengths. This is donebecause the solutions themselves contain nodal points from interfering waves, where theamplitude is very small, and are hence not directly suitable as reference value. Very smallrelative errors are obtained (except directly at the source point), ranging from less than0.01 over most of the domain to 0.05 or 0.08 at isolated spots where the absolute amplitudeis small. 7. Application in multigrid based solvers
The last few years there have been several interesting developments in multigrid methodsfor Helmholtz equations. Different two-grid methods with inexact coarse level solvershave been studied in [5] and [25]. In [5] a number of iterations of shifted Laplacianpreconditioned Krylov solver [11] is used as coarse level solver. The method of [25] is basedon the multigrid method in [26] with a double sweep domain decomposition preconditioner[24] as coarse level solver. The multigrid method with exact coarse level solver was studied
DISPERSION MINIMIZING SCHEME FOR THE 3-D HELMHOLTZ EQUATION 27 solution at 50 Hzx(m) y ( m ) Figure 8.
Solution from a point source at 50 Hz(a) (b) x(m) y ( m ) reference amplitude at 50 Hz x(m) y ( m ) reference amplitude at 100 Hz (c) (d) x(m) y ( m ) relative error at 50 Hz x(m) y ( m ) relative error at 100 Hz Figure 9.
Reference values for 50 Hz (a) and for 100 Hz (b). Relativeerrors in the solutions for 50 Hz (c) and for 100 Hz (d).in [26]. There it was shown that the convergence can be strongly improved when phaseslowness differences between the fine and coarse scale operators are minimized. For thispurpose, optimized finite differences were used at the coarse level, and good convergencewas obtained for meshes with downto three points per wavelength at the coarse level. Forstandard choices of the coarse level discretization it was found that about 10 points perwavelength at the coarse level were needed to have good convergence.In [26] standard second order finite differences were used as the fine level. Becauseof the relatively large phase slowness errors of this method, the coarse level optimizedfinite difference method had to be constructed specifically to match the phase slownesssof second order finite differences, instead of matching the true phase slowness. A betterchoice is to use method with small phase slowness errors at the fine level and at the coarselevels. Here we will use IOFD at all levels. These experiments do not involve the operator (a)
Marmousi modelx(m) y ( m ) (b) − D slice of the SEG − EAGE salt modelx(m) y ( m ) Figure 10.
Velocity models for the 2-D two-grid experiments. (a) Mar-mousi (b) 2-D slice of the SEG-EAGE salt model. Q . The operator P is used directly as coarse level discretization and at the fine level weare only interested in solving equation (111).In the first set of computational results of this section we will show that this resultsin good convergence of the multigrid method with exact coarse level solver. In a secondexample we will study the multigrid method with inexact coarse level solver of [25].In our study of the convergence when using the exact coarse level solver we are againinterested in examples with wave propagation over hundreds of wavelengths. Therefore,these experiments are done in two dimensions. For background on multigrid methods,see [28]. As in [26] most of the components of the multigrid method are standard. Fullweighting restriction and prolongation operators are used. As smoother, an ω -Jacobimethod is used. We found that ω = 0 . ν = 4 (the number of pre- and postsmoothingparameters) are good choices of parameters. In these experiments we used a conventionalabsorbing boundary layer to simulate an unbounded domain.We studied the convergence as a function of the number of points per wavelength forthree velocity models: A constant model, the Marmousi model and a slice of the 3-D SEG-EAGE salt model. The latter two models are displayed in Figure 10. The parameters ofthe examples and the observed number of iterations to reduce the residual by 10 are givenin 5. At downto three points per wavelength the method behaved well. At 2.5 ppw coarselevel the method still converged, but the number of iterations increased substantially, andalso became more sensitive to the problem size (which was apparent from smaller scaleexperiments not included in the table).Note that the application in multigrid methods is quite different from the applicationas fine level discretization. The method is used at coarser meshes (at three points perwavelength the direct application will in general lead to too large errors). Also, multi-grid solvers using IOFD at coarse levels may be developed for other types of fine leveldiscretizations, as long as they use a regular mesh.We now turn to a multigrid method with an inexact coarse level solver. Such methodsare used because in three dimensions it is often too expensive to compute the exact solu-tion. These methods are currently some of the fastest solvers for large problems that arein the literature [5, 25].Because we are interested in coarse meshes, such as six points per wavelength based onthe previous examples, it is a priori not clear that the above mentioned solvers perform DISPERSION MINIMIZING SCHEME FOR THE 3-D HELMHOLTZ EQUATION 29 constant Marmousi salt model2400 × ×
750 2700 × Table 5.
Iterations required for a two-grid method using IOFD discretiza-tion at the fine and coarse level as a function of the number of points perwavelength (ppw)well. Like many solvers in the literature, they were tested for problems with at leastten mesh points per wavelength. They cannot be assumed to converge as well for largerfrequencies, because multigrid convergence depends on frequency, and the same is true forthe shifted Laplacian preconditioner [9]. For the double sweep domain decomposition it isunclear how the frequency affects the convergence, but a priori it also cannot be assumedto be independent of the frequency.This raises the question whether we can actually obtain a gain in efficiency by goingto coarser meshes. The purpose of the next example is to show that this indeed the case,and to generally show that IOFD can perform well with the solver of [25].In the following example we will test the method of [25], which is a two-grid methodusing an inexact coarse level solver given by a double sweep domain decomposition pre-conditioner (see [24]). The method is modified to use IOFD at both the fine and the coarselevels of the two-grid method. We will take the SEG-EAGE Salt Model as an example,similarly as in [25]. In addition to changing the discretization method we will increase thefrequency by a factor , so that a minimum of six points per wavelength is used, a regimewhich has not been tested before for this method. If convergence and cost per degree offreedom would stay constant, there would be an improvement in the cost by a factor ofover (cid:0) (cid:1) ≈ .
62 (more than this because cost grows somewhat faster than linear withproblem size).The original SEG-EAGE salt model is of size 13500 x 13500 x 4200 meter, discretizedwith 20 m grid spacing. We apply the method just described to solve the Helmholtzequation with this velocity model and random or point sources as right hand sides at fourdifferent frequencies from 6 .
25 to 12 . ν = 3 for the number of pre- and postsmoothingsteps and ω = 0 .
65 in the ω -Jacobi method. Computations were done on the Lisa clusterat Surfsara, described already in section 6, using the implementation described in [25]. Amaximum of 16 nodes were used in parallel for these computations.The algorithm is set up to solve for multiple right hand sides simultaneously. In thetable of results, the computation time per right hand side is given. In Table 6 someparameters are given, together with the computation time and iteration count to reducethe residual by 10 − . As illustration, plots of a solution are given in Figure 12. It can beobserved that the cost increases very little compared to the results of [25], even thoughfrequencies are increased by a factor 5 /
3. Some increase in cost can be expected, becausethe discrete Helmholtz operator using second order finite differences is cheaper to applythan the one using a compact 27-point stencil. Hence reducing the number of points perwavelength in the mesh can indeed lead to corresponding savings in computation time. (a) (b) x z SEG EAGE salt model 0 2000 4000 6000 8000 10000 1200005001000150020002500300035004000 150020002500300035004000 y z SEG EAGE salt model 0 2000 4000 6000 8000 10000 1200005001000150020002500300035004000 150020002500300035004000
Figure 11.
SEG-EAGE salt velocity model: (a) ( x, z ) slice at y = 6740m (b) ( y, z ) slice at x = 6740 m.(a) (b) x z Wavefield in SEG EAGE salt model0 2000 4000 6000 8000 10000 1200005001000150020002500300035004000 y z Wavefield in SEG EAGE salt model0 2000 4000 6000 8000 10000 1200005001000150020002500300035004000
Figure 12.
Solution to the Helmholtz equation at 12 . x, z ) sliceat y = 6740 m (b) ( y, z ) slice at x = 6740 m.frequency 6.25 7.87 9.91 12.5size 338x338x106 426x426x132 536x536x166 676x676x210 . · . · . · . · cores 32 64 128 256 Table 6.
Computation times and iteration counts for the SEG-EAGE SaltModel example.As mentioned, the methods of [5] and [25] are some of the fastest currently in theliterature. Comparing with these results we see a significant improvement. For examplein [5] the SEG-EAGE salt model problem was solved at 10 Hz in 270 seconds on 256 coresof an IBM BG/P machine (with the residual reduced by a factor 10 instead of 10 in ourcase). Here we solve the problem at 9.91 Hz using 128 cores in 45 seconds per right handsides (179 seconds for four right hand sides), a clear improvement .8. Discussion
Here we summarize some of the conclusions and further discuss the results.Using the results presented one can make a case for the use of coarse meshes using aminimum of five or six points per wavelength in time harmonic wave simulations in case k is smooth. This idea is not new, in the exploration geophysics community it appears tobe quite common. However, we found that the methods that have been proposed for thispurpose in [16] and [21] can be expected to give substantial phase errors in simulationsof large distance wave propagation. By using the new IOFD method (in two or three On the other hand the method of [5] uses less memory and has been applied to larger examples thanwe have shown here. Furthermore no full comparison including accuracy was made.
DISPERSION MINIMIZING SCHEME FOR THE 3-D HELMHOLTZ EQUATION 31 dimensions) or the QS-FEM method (in two dimensions only), phase errors can be mademuch smaller.In k ( x ) has strong gradients, one can expect that, at least locally where ∇ k is large, finermeshes and/or different discretizations are needed to obtain accurate solutions. Large gra-dients lead to reflections. Typically finer meshes are needed to model these accurately. Forone reason this is because finer discretizations of k are needed, since linearized scatteringtheory shows that reflected waves are associated with perturbations in the medium veloc-ity with wave vectors of length up to 2 k (where here k refers to the background velocityaround which the linearization is applied). However, in this case the multigrid approachdiscussed in section 7 can still be useful. It has been applied in successfully in exampleswith discontinuities. This suggests to do further research on multigrid approaches withcompact finite difference method at the coarse level and other discretizations at the finelevel, including methods with local refinement. A similar argument can be held for thediscretization of the right hand side f in the equation ( − ∆ − k ( x ) ) u = f . For rapidlyvarying functions f , finer meshes may be needed at least locally where the rapid variationsoccur.When applied in inversion algorithms IOFD and QS-FEM are somewhat more com-plicated than the methods of [16] and [21], because the operatore depends in a morecomplicated fashion on the coefficients, which means it is more complicated to computethe derivative of the finite difference operator with respect to the medium coefficients.Due to the use of Hermitian interpolation, these derivatives are however continuous forour IOFD method. References [1] S. Alinhac and P. G´erard.
Pseudo-differential operators and the Nash-Moser theorem , volume 82 of
Graduate Studies in Mathematics . American Mathematical Society, Providence, RI, 2007. Translatedfrom the 1991 French original by Stephen S. Wilson.[2] P. R. Amestoy, I. S. Duff, J. Koster, and J.-Y. L’Excellent. A fully asynchronous multifrontal solverusing distributed dynamic scheduling.
SIAM Journal on Matrix Analysis and Applications , 23(1):15–41, 2001.[3] I. Babuˇska, F. Ihlenburg, E. T. Paik, and S. A. Sauter. A generalized finite element method forsolving the Helmholtz equation in two dimensions with minimal pollution.
Comput. Methods Appl.Mech. Engrg. , 128(3-4):325–359, 1995.[4] J.-P. Berenger. A perfectly matched layer for the absorption of electromagnetic waves.
Journal ofComputational Physics , 114(2):185–200, 1994.[5] H. Calandra, S. Gratton, X. Pinel, and X. Vasseur. An improved two-grid preconditioner for thesolution of three-dimensional Helmholtz problems in heterogeneous media.
Numer. Linear AlgebraAppl. , 20(4):663–688, 2013.[6] Z. Chen, D. Cheng, and T. Wu. A dispersion minimizing finite difference scheme and preconditionedsolver for the 3D Helmholtz equation.
Journal of Computational Physics , 231(24):8152 – 8175, 2012.[7] W. C. Chew and W. H. Weedon. A 3D perfectly matched medium from modified Maxwell’s equationswith stretched coordinates.
Microwave and optical technology letters , 7(13):599–604, 1994.[8] G. C. Cohen.
Higher-order numerical methods for transient wave equations . Scientific Computation.Springer-Verlag, Berlin, 2002. With a foreword by R. Glowinski.[9] S. Cools and W. Vanroose. Local Fourier analysis of the complex shifted Laplacian preconditioner forHelmholtz problems.
Numer. Linear Algebra Appl. , 20(4):575–597, 2013.[10] J. J. Duistermaat.
Fourier integral operators , volume 130 of
Progress in Mathematics . Birkh¨auserBoston, Inc., Boston, MA, 1996.[11] Y. A. Erlangga, C. W. Oosterlee, and C. Vuik. A novel multigrid based preconditioner for heteroge-neous Helmholtz problems.
SIAM J. Sci. Comput. , 27(4):1471–1492 (electronic), 2006.[12] G. Grubb.
Distributions and operators , volume 252 of
Graduate Texts in Mathematics . Springer, NewYork, 2009.[13] I. Harari and E. Turkel. Accurate finite difference methods for time-harmonic wave propagation.
J.Comput. Phys. , 119(2):252–270, 1995. [14] L. H¨ormander.
The analysis of linear partial differential operators. III , volume 274 of
Grundlehren derMathematischen Wissenschaften [Fundamental Principles of Mathematical Sciences] . Springer-Verlag,Berlin, 1985. Pseudodifferential operators.[15] F. Ihlenburg.
Finite element analysis of acoustic scattering , volume 132 of
Applied MathematicalSciences . Springer-Verlag, New York, 1998.[16] Jo, Churl-Hyun and Shin, Changsoo and Suh, Jung Hee. An optimal 9-point, finite-difference,frequency-space, 2-D scalar wave extrapolator.
Geophysics , 61(2):529–537, 1996.[17] S. G. Johnson. Notes on perfectly matched layers. http://math.mit.edu/ stevenj/18.369/pml.pdf, 2010.[18] H. Knibbe, W. A. Mulder, C. W. Oosterlee, and C. Vuik. Closing the performance gap between aniterative frequency-domain solver and an explicit time-domain scheme for 3D migration on parallelarchitectures.
GEOPHYSICS , 79(2):S47–S61, MAR-APR 2014.[19] M. J. Lighthill. Studies on magneto-hydrodynamic waves and other anisotropic wave motions.
Philos.Trans. Roy. Soc. London Ser. A , 252:397–430, 1960.[20] J. M. Melenk and S. Sauter. Wavenumber explicit convergence analysis for galerkin discretizations ofthe helmholtz equation.
SIAM J. Numer. Anal. , 49(3):1210–1243, June 2011.[21] S. Operto, J. Virieux, P. Amestoy, J.-Y. L’Excellent, L. Giraud, and H. B. H. Ali. 3D finite-differencefrequency-domain modeling of visco-acoustic wave propagation using a massively parallel direct solver:A feasibility study.
Geophysics , 72(5, S):SM195–SM211, 2007.[22] J. Poulson, B. Engquist, S. Li, and L. Ying. A parallel sweeping preconditioner for heterogeneous 3DHelmholtz equations.
SIAM J. Sci. Comput. , 35(3):C194–C212, 2013.[23] C. D. Riyanti, A. Kononov, Y. A. Erlangga, C. Vuik, C. W. Oosterlee, R.-E. Plessix, and W. A. Mulder.A parallel multigrid-based preconditioner for the 3d heterogeneous high-frequency helmholtz equation.
Journal of Computational Physics , 224(1):431 – 448, 2007. Special Issue Dedicated to Professor PietWesseling on the occasion of his retirement from Delft University of Technology.[24] C. C. Stolk. A rapidly converging domain decomposition method for the Helmholtz equation.
J.Comput. Phys. , 241:240–252, 2013.[25] C. C. Stolk. A two-grid accelerated sweeping preconditioner for the Helmholtz equation, 2014.Arxiv.org/abs/1412.0464.[26] C. C. Stolk, M. Ahmed, and S. K. Bhowmik. A multigrid method for the helmholtz equation withoptimized coarse grid corrections.
SIAM Journal on Scientific Computing , 36(6):A2819–A2841, 2014.[27] G. Sutmann. Compact finite difference schemes of sixth order for the Helmholtz equation.
J. Comput.Appl. Math. , 203(1):15–31, 2007.[28] U. Trottenberg, C. W. Oosterlee, and A. Sch¨uller.
Multigrid . Academic Press Inc., San Diego, CA,2001. With contributions by A. Brandt, P. Oswald and K. St¨uben.[29] E. Turkel, D. Gordon, R. Gordon, and S. Tsynkov. Compact 2D and 3D sixth order schemes for theHelmholtz equation with variable wave number.
Journal of Computational Physics , 232(1):272 – 287,2013.
Appendix A. Phase slowness computations for finite element methods
In a periodically repeating setting, which is the case for finite elements with N ≥ p, q, r divisible by N . For such operators we consider the Bloch waves(121) u l,m,n = e iξ · x l,m,n v l,m,n where v l,m,n is periodic with shifts ( pN, qN, rN ), p, q, r integers. For given ξ , the action ofan operator A with these symmetries is given by an N × N matrix acting on the v l,m,n (in three dimensions) for 0 ≤ l, m, n ≤ N −
1. We need to find the vectors ξ for whichthere is a zero eigenvalue. However not all zero eigenvectors correspond to plane waveswith wave vectors ξ , because of the presence of v l,m,n . In general v l,m,n can correspond toa linear combination of plane waves with wave vectors given by ( p πNh , q πNh , r πNh ), where p, q, r are integers. Assuming that some ξ corresponds to a simple zero eigenvalue andthat u l,m,n is close to a plane wave (which is often the case because the eigenfunctionsof the continuous operator are plane waves and the operator A is a good approximationof the continuous operator), these integers p, q, r can be determined modulo N , and thewave vector associated with an element of the zero set of A can be determined.In the computations we will take a somewhat different approach. We will computeall eigenvalues, and then only consider the one whose eigenvector v l,m,n is most closely DISPERSION MINIMIZING SCHEME FOR THE 3-D HELMHOLTZ EQUATION 33 correlated with (has the in absolute value largest inner product with) the constant function˜ v l,m,n = 1. We will say we have found a phase velocity vector at some ξ if this eigenvalueis zero. This approach has some limitations, but a more extensive study of this topic fallsoutside the scope of this paper. For standard finite elements and k, hk, h