1 Introduction
The Helmholtz equation is used to model the propagation of a wave within a heterogeneous medium. Its acoustic version is given by
(1) 
where
is the Fourier transform of the wave’s
pressure field, is the angular frequency, is the “slowness” of the wave in the medium (the inverse of the wave velocity), and is the density of the medium. The righthandside incorporates sources into the equation. The equation is discretized on a finite domain and is accompanied with some absorbing boundary conditions that mimic the propagation of a wave in an open domain. This is usually achieved by some complex valued absorbing boundary layer [16, 31, 46], which is related to modeling the attenuation of the wavefield.The acoustic equation (1) is usually discretized by a finitedifference scheme on a regular grid, resulting in a large and indefinite linear system, which is complexvalued due to the absorbing boundary conditions and possible attenuation. If the frequency (or the wavenumber ) is high, the problem requires a very fine mesh and a large number of unknowns [3, 25]. In this common case, solving the discretized equation at large scale 3D scenarios is challenging, and is still considered to be an open problem.
One of the main applications that include the Helmholtz equation is Full Waveform Inversion [39, 60, 34, 50]
, which is a process used to estimate the wave velocity and rock structure of the earth’s subsurface. The inversion process (in the frequency domain) includes many repeated solutions of Helmholtz equations for modeling the wave propagation. These solutions are used to iteratively estimate the unknown wave velocity in the earth’s subsurface. However, because the earth is an elastic medium, the acoustic equation in (
1) does not fully capture the physics of the wave propagation, and research is advancing towards elastic waveform inversion, in which the elastic Helmholtz equation is solved for modeling the wave propagation [18, 9, 8, 5, 27].The elastic Helmholtz equation, which we formulate later, is a system of partial differential equations (PDEs). While (
1) models pressure waves only, the elastic Helmholtz equation also models shear waves. Similarly to (1), the linear system that results from discretizing the elastic equation is indefinite and complex valued. Moreover, because the equation is a system of PDEs, then the associated linear system is three times larger (in 3D) than the acoustic one for the same mesh size. In addition, the discretization requires the mesh to be finer than in the acoustic (1), because the modeled shear waves have higher wavenumber than the pressure waves [33]. Altogether, we get a huge linear system which is more difficult to solve than the acoustic one and an iterative method is required for its solution. However, while the solution of the acoustic equation (1) has been heavily studied in the literature with a variety of methods [17, 38, 25, 32, 6, 35, 23, 51, 49, 21, 11], the elastic version has very few available iterative solvers known to us. One recent “elastic solver” is [30], which is an extension of [23] to the elastic case. This method, which involves a hybrid parallel Kaczmarz preconditioner, is quite simple and generic, and hence requires a lot of iterations to solve the system at large scales.One of the most common solvers for the discretized acoustic equation (1) is the shifted Laplacian multigrid method [20, 19, 55, 36, 53, 10, 1, 14, 13, 40], where an attenuated version of (1) is used as a preconditioner for the true system inside a Krylov method. The attenuated system, which is the same system with a complex shift, can be easily solved by multigrid, if the attenuation is high enough. However, as we add more attenuation, the efficiency of preconditioner deteriorates. This is a tradeoff that methods try to balance.
The shifted Laplacian multigrid approach seems to be naturally extendable to the elastic case. However, experiments show that the standard shifted Laplacian multigrid is not efficient for the elastic Helmholtz equation, and some specialized treatment is necessary. Indeed, the recent [41] suggests such a multigrid method, using linerelaxation instead of pointwise relaxation. However, this method is only presented for 2D problems and does not seem to achieve the same efficiency compared to the acoustic case. In particular, the authors use a significantly higher attenuation (shift) parameter than what is usually used in the acoustic case, and the line relaxations in 2D extend to relatively expensive plane relaxations in 3D. These two properties leave room for improvement.
In this work, we present a new shifted Laplacian multigrid method for the elastic Helmholtz equation. Our method adopts approaches that are suitable for linear elasticity to better treat the elliptic “elastic part” of the elastic Helmholtz equation, which in some sense plays the same role as the weighted Laplacian in (1). This elliptic part is essentially the elasticity operator, and is known to cause difficulties to standard multigrid methods in cases of nearly incompressible material. The shifted Laplacian method is no exception, and to solve (2) using multigrid, we require to apply the mechanisms for both the elasticity and the acoustic Helmholtz equation together. To this end, we write the elastic equation using a mixed formulation [22] and use a local cellwise “Vanka” relaxation to treat the elastic part of the elastic Helmholtz equation [61]. The indefiniteness of the problem is treated by shifted Laplacian in the same way that the indefiniteness of (1) is treated. We demonstrate that our method performs similarly to standard shifted Laplacian for (1), only with respect to the shear wavenumber instead of the pressure wavenumber.
Our paper is organized as follows: in the next two subsections we present the elastic Helmholtz equation, discuss its discretization, and briefly present the general shifted Laplacian method. In Section 2 we present our method, and discuss the relation between the acoustic and elastic equations in the case of fully incompressible material. Finally in Section 3 we demonstrate the properties and efficiency of our method in a few numerical examples in two and three dimensions.
1.1 Problem formulation and discretization
The elastic Helmholtz equation has several formulations. Here we focus on the equation in isotropic medium, which is formulated by either of the following equivalent equations
(2)  
The unknown
is a displacement vector which has three components at each location in the domain.
and are the Lame parameters, and is the density of the medium as in (1). These parameters determine the pressure and shear wave velocities by , and , respectively [29]. The term(3) 
is the strain tensor (factored by two) and is a symmetric tensor of nine scalar quantities in 3D. The divergence that operates on the strain tensor (
3) in (2) essentially operates on each row or column of the tensor, making each row/column into a scalar quantity. The term is the weighted diffusion operator applied on each of the components of the vector separately. In the case of , (2) become the elasticity operator. In the case where the material is incompressible. Then, as we show later, the elastic equation can be reduced to the acoustic equation (1), modeling only pressure waves, with as the pressure wave velocity. Other richer elastic formulations may include more parameters than and , e.g. the orthorhombic formulation has 9 parameters and also models anisotropy [30]. In principle, the method that we present in this paper is suitable for anisotropic cases as well as long as the anisotropy is not too strong. If the anisotropy is strong, our method will require adaptations similar to those that are needed for the standard shifted Laplacian method for an anisotropic Laplacian operator in (1), e.g. semicoarsening. Such extensions are beyond the scope of this paper.For simplicity we focus on a discretization using finitedifferences schemes on a regular rectangular mesh. To discretize (2) one has to choose between two main approaches:

A nodebased discretization where all the components and are placed at the same locations on the nodes of the grid cell.

A staggered grid approach where the components and are placed on different faces or edges of the grid cell.
The advantage of the nodal approach (option 1) is that one can formulate highorder or optimally weighted secondorder discretizations using a compact 27point stencil [26, 48, 24] at each component. The reasons for using this approach are (1) such discretizations leads to low dispersion error and (2) the 27point stencil is a favorable sparsity pattern for numerical solvers. This pattern leads to fast memory access to a vector when applying a matrixvector multiplication, and low fillin when using direct solvers such as MUMPS [2] or PARDISO [44]. For the same reasons, similar nodal discretization schemes are also used for discretizing the acoustic (1) in [37, 45, 47, 54]. On the other hand, the disadvantage of the nodal approach for discretizing the elastic equation (2) is that there are relatively large errors in cases of nearly incompressible materials when [41]. The situation in the staggered discretization (option 2) is the opposite from the one in the nodal discretization. On the one hand, the staggered discretization is stable for nearly incompressible materials [59, 28]. On the other hand, compact highorder finite difference discretizations are currently not available. For example, the fourth order schemes of [28, 29] create stencils that are wider than the 27point stencil block, because the gradient itself is approximated by a fourth order operator. This leads to relatively high memory access time in matrixvector products, and relatively high fillin using direct methods. Still, the high order schemes are important for minimizing the dispersion error that is significant in standard second order schemes, if the mesh is not fine enough.
In this work we focus on a multigrid solver for (2), and we use a standard staggered grid discretization, illustrated in Figure 1. This discretization is similar to the one used in [59, 41], and suggested in [52] to discretize similar systems of PDEs. We place the components of on the faces of the cell, and denote its discrete components in boldface, i.e., the vector has the discrete values of on the mesh. Similarly, are placed on the cell center and their discrete vectors are denoted by .
We discretize the first derivatives using central differences
This way, the numerical derivative of each component with respect to each variable ends up at a different place in the cell. For example, the first derivative of with respect to ends up in the cell center, at the same location as the discrete . The derivative of with respect to ends up on the middle of an edge, and on this same edge we also get the first derivative of with respect to . This property is necessary, because the discrete components of the strain tensor (3) need to end up at the same locations.
The linear system that results from discretizing the bottom formulation in (2) is given by
(4) 
where is the cellcentered gradient operator that operates from the cell centers to the faces, where the components are placed. The matrix is used for the discrete divergence operator. The operator generates a diagonal matrix with the values of on the cell centers. The operator is a block diagonal gradient matrix, which includes three gradient matrices on its diagonal—one for each of the components —using central difference schemes. This operator computes one of the terms in Eq. (3) (which one depends only on the ordering of the strain tensor). averages the values of the cellcentered to the edges, and creates a diagonal matrix to hold the averaged values. Altogether, the term ends up being the weighted Laplacian operator which is applied on each of the components separately. We choose the mass matrix
(5) 
where is a diagonal matrix with the averaged values of a cellcentered onto the cell faces. The symbol is the Hadamard product, and the vector is a physical attenuation vector. We also use to incorporate the absorbing boundary conditions, using a function that quadratically goes from zero to one towards the domain boundaries [12]. The attenuation can be equivalently modeled by using complex frequencydependent Lame parameters [48], which can also be used to model different attenuation factors for the shear and pressure wave velocities [24]. We place on the cellcenter and assume that the physical attenuation is very small.
One nice feature of the staggered discretization is that both the equations in (2) are equal even after the discretization (up to boundary conditions). That is, we get two identical linear systems for each one of the formulations, up to small differences that result from averaging the coefficients and boundary conditions only. That is not the case with nodal discretizations. Finally, it is worthy to note that we use this second order scheme, as in [59, 41], to demonstrate our solver, but expect that it will also be effective for the highorder staggered discretizations [28, 29]. That is similarly to [20] and [55] which show that the shifted Laplacian methods works with similar efficiency for both second order and highorder nodal schemes for the acoustic Helmholtz case.
1.2 Multigrid methods and the shifted Laplacian framework
In this section we describe the general shifted Laplacian multigrid framework that we adopt in this paper. Multigrid methods aim at solving linear systems
(6) 
iteratively by using two complementary processes. The first process is the relaxation, which is obtained by a standard local iterative method like Jacobi or GaussSeidel. Such methods are typically effective at reducing only part of the error in the iterative solution process. The other part of the error, called “algebraically smooth”, is not reduced well by the relaxation and is typically defined by vectors s.t.
(7) 
To reduce such errors, multigrid methods use a “coarse grid correction”. In this correction, the error for some iterate is estimated by solving a coarser system
where the matrix is an approximation of the matrix on a coarser mesh (the subscript denotes coarse components). The matrix
is the so called prolongation operator that is used to restrict the residual onto the coarser grid, and to interpolate the solution of the coarse system,
, back to the fine grid:(8) 
The coarse operator can be obtained by either rediscretizing the problem on a coarser grid or by the Galerkin operator
(9) 
Algorithm 1 summarizes the process using two grids. By treating the coarse problem recursively, we obtain the multigrid Vcycle, and by treating the coarse problem recursively twice (by two recursive calls) we obtain a Wcycle. This multigrid process is effective if all smooth errors satisfying (7) are represented in the range of the prolongation . For more information see [7, 52] and references therein.
To solve Helmholtz problems like (1) efficiently, the process above requires some modification, like in the shifted Laplacian framework. To apply this framework, one introduces a shifted matrix
(10) 
where is a matrix defined by a discretization of a Helmholtz operator, is some mass matrix, and is a shifting parameter. Usually, is defined as a mass matrix that is used for modeling attenuation in Helmholtz systems. In the acoustic case , and in our elastic case . The advantage of the shifted system is that it can be efficiently solved by multigrid methods. Hence, in the shifted Laplacian framework, the shifted Helmholtz matrix (10) is used as a preconditioner for a Helmholtz linear system (6) inside a suitable Krylov method like (flexible) GMRES [43] or BiCGSTAB [56]. The preconditioning is obtained by approximately inverting the shifted matrix (10) using a multigrid cycle.
When modeling waves with an attenuated Helmholtz problem (with instead of ), the waves decay rapidly if is large. As we add more attenuation (larger in (10)), we can invert the shifted matrix more easily, but the performance of the shifted Laplacian preconditioner deteriorates. This is a tradeoff that methods try to balance, and the common compromise chosen in [20] for the acoustic equation is to use . [10] and [50] suggest using less attenuation, but invest more effort in the multigrid cycles. In [41] the elastic equation (4) is solved, and the authors use a high shift parameter () in (10), which results in a much less efficient preconditioner regardless of the multigrid ability to invert the shifted matrix.
The prolongation is usually chosen to be a bilinear interpolation operator, or an operator induced prolongation like AMG [20]. As relaxation, the damped Jacobi method [50] or the GMRES method [15, 10, 14], are often chosen for the acoustic case. For the elastic case, linerelaxation may be used [41]. Beside the relaxation, one has to choose the number of levels used in the multigrid hierarchy. Unlike other scenarios, the algebraically smooth error modes of the Helmholtz operator are not represented well on a very coarse grid because of their signchanging patterns. As a result, the performance of the solver deteriorates as we use more levels at high frequency. The common choice is using three levels only [10, 50, 19], but invest more work on solving the second grid. However, using only three levels is problematic in large 3D cases, as we get a rather large coarsest grid problem, which is usually solved by a direct solver. Factorizing these matrices is highly memory consuming and limiting, and it is even more so for the elastic systems. The work in [10] suggests an inexact solution of the coarsest grid using GMRES (for the acoustic case). We elaborate on our choices of relaxation, interpolation and other parameters when we present our method for solving (4) in the next section.
2 Shifted Laplacian for the elastic Helmholtz equation
It is clear that if we wish to solve a Helmholtz equation, whether (1) or (2), we first need to be able to solve the equation for a low frequency . In (1), we are left with the elliptic weighted Poisson equation, which is considered to be the “bread and butter” of multigrid methods and can be easily solved even in cases of jumping coefficients and anisotropy. On the other hand, if we set in (2) we get the linear elasticity equation which is also elliptic, but more difficult to solve than the Poisson equation. Thus, the idea that guides us is: if we wish to solve (2) or invert its shifted version, then our solver must be able to handle the elasticity problem efficiently.
Linear elasticity problems can be solved by multigrid methods, but they require some special treatment in the nearly incompressible case when . In this case, the linear elasticity problem has a dominating graddiv operator, which has a rich nullspace. Indeed, because , then any vector function which is a curl of another vector function is in the nullspace of the graddiv operator [22]
This equality also holds in discrete space using staggered discretization. This requires a special treatment, because a simple prolongation operator cannot not approximate this rich null space well in its range, causing simple multigrid methods to be inefficient.
There are several approaches to handle this rich near nullspace. One is by using a prolongation that has this nullspace in range, e.g. a prolongation based on smoothed aggregation which includes all the rigid body modes as basis functions [57]. This results in a rather large coarse grid matrix, but does not require further modifications to the relaxation or other multigrid ingredients. A different family of approaches involves rather standard transfer operators, but also involves reformulating the system (4) into an equivalent one, called “mixed formulation” [22, 62]. The mixed formulation, which is the approach that we choose in this work, is achieved by introducing a new pressure variable . In discrete form
(11) 
where the operator is the cellcentered diagonal matrix operator defined right after (4). We then reformulate the linear system (4) as the coupled system
(12) 
This linear system is equivalent to (4), because by construction, the system (4) is the Schur complement of (12) when eliminating the block.
Using the reformulated system is the first step in tackling the problem. The second step is to use a special relaxation scheme. Here there are two options: either the socalled “distributive relaxation” [22, 62], or the cellwise “Vanka” relaxation [61]. The latter was originally developed for the Stokes equation [58], which has a similar structure as (12) but with a zero block for the variable. Both options are suitable for handling the problem, and in this work we choose the cellwise relaxation [61] which in our opinion is simpler to implement and to parallelize in multicore computations than the distributive relaxation.
In the cellwise relaxation, we sweep through all the cells in the domain, and for each cell we invert the local matrix composed of the block of (12) for that cell. In 3D, each cell has two variables for each displacement component on each of the faces, and the pressure variable at cellcenter—a total of 7 variables, that leads to a block of the matrix. To apply this in parallel, we use the redblack cell wise relaxation, where the all the cells in the domain are divided according to a checkerboard pattern into “reds” and “blacks”, and we first simultaneously sweep through all the red cells, and then simultaneously sweep through all the black cells. Figure 2 shows an example of the cellwise relaxation in two dimensions.
Summary of the method
To solve the elastic Helmholtz equation we first reformulate the system (4) to an equivalent mixed formulation system (12). In our multigrid cycles we use standard bilinear transfer operators that are suitable for facebased staggered discretization [52]. As relaxation we use redblack cellwise relaxation with a damping parameter of 0.5. The multigrid hierarchy is defined for a shifted version of (12
), which is added with the zeropadded shift matrix
That is, the artificial attenuation that we add involves only the block, just as we would have done for the original formulation (4). The attenuation is not involved with the pressure variable resulting from the mixed formulation. Finally, the system (12) is solved using a Krylov method with a multigrid cycle for the shifted operator as a preconditioner.
The relation between the acoustic and elastic Helmholtz equations
As noted before, if , then the elastic equation models only pressure waves. Interestingly, the formulation (12) encapsulates this. The pressure variable introduced in (11) is in fact the same pressure waveform function that appears in the acoustic (1), up to a diagonal scaling. Indeed, if we set and ignore the attenuation parameter in (5), then the Schur complement of (12) when eliminating the block ends up as a cellcentered discretization of (1). This means that if we have a problem that is part acoustic and part elastic we can formulate the problem using one discrete system (12), but treat both parts separately using the conversions in (11) and
that is coming from the top block in (12). This will result in a domaindecomposition type of solver, where the acoustic part is solved using standard shifted Laplacian, and the elastic part with elastic shifted Laplacian. The acoustic solver is obviously cheaper than the elastic version. Such scenarios of mixed elastic and acoustic media are common in marine full waveform inversion, where the sea is and acoustic medium and the rock under it is an elastic medium.
A remark on implementation:
The cellwise relaxation involves sweeping through all cells, and inverting an essentially dense matrix for each one. That may be a costly operation. Instead, in our implementation we extract these local matrices and invert them at the setup phase. In the solve phase which includes quite a few applications of the cellwise relaxation, we only multiply the values of the inverted matrices instead of extracting and inverting the local matrices onthefly. This results in a storage of variables in 3D (note that the number of cells in approximately one fourth of the matrix variables). To handle the storage issue, we convert and save the inverted matrices in a low 16bit half precision. This results in quite fast application of the cellwise relaxation in the price of a moderate storage requirement. Other ingredients like the coarsest grid factorization, and the vectors needed for the Krylov method are much more memory consuming in our experience.
3 Numerical Results
In this section we demonstrate the shifted Laplacian solver for the elastic Helmholtz equation. We present examples that appear in geophysical applications, where typically the length of the domain is quite high (about 20km) compared to its depth (about 5km). We first perform a few experimental comparisons in 2D, and then perform 3D experiments. Note that the 2D experiments are done only for illustrating the behavior of the solver for 3D problems at those scales. The first experiment compares the shifted Laplacian method based on standard formulation vs the mixed formulation. The second one compares the performance of shifted Laplacian for the acoustic equation using the shear velocity versus shifted Laplacian for the elastic equation using mixed formulation. The third one shows the effect using three vs. four levels in the multigrid cycles. The last experiment shows some 3D performance. Our code is written in the Julia language [4], and is written as part of the jInv.jl package [42]. Using this package, our code can be easily used as a forward solver for three dimensional elastic full waveform inversion in the frequency domain.
In all the examples we use the preconditioned GMRES(5) Krylov solver [43], and seek a solution with relative residual accuracy of , starting from a zero initial guess. The right hand side is chosen to be a point source located at middle of the top row of the domain, similarly to Fig. 3. Unless stated otherwise, we use the shifted operator (10) as preconditioner, using a shift parameter of 0.2. This shift parameter was used in [50, 51] for solving the acoustic equation (1). In all the experiments that concern the solution of the reformulated system (12), we use 3level cycles using redblack cellwise relaxation, with a damping parameter of 0.5. In the multigrid framework, we define our coarse grid problems by the Galerkin product (9). These are our ”default” parameters that we found to be the most effective. In all the examples, we use an absorbing boundary layer of 20 cells, and add a small artificial attenuation by adding to in (5).
3.1 Standard vs. mixed formulations
Shifted Laplacian: standard vs. mixed formulation  

= 1  = 2  = 4  = 8  = 16  
= 0.25  = 0.33  = 0.4  = 0.44  = 0.47  
35 (35)  35 (35)  35 (35)  38 (35)  125 (35)  
41 (38)  43 (38)  45 (37)  51 (37)  122 (37)  
58 (44)  66 (44)  79 (44)  110 (44)  (43 )  
110 (67)  138 (69)  179 (69)  (71)  (71)  
220 (108)  (112)  (117)  (115)  (120)  
In our first experiment we use a twodimensional constant coefficients example, and monitor the influence of and on the performance of the shifted Laplacian method, while keeping and fixed. The case of nearly incompressible material correspond to a case where , or when the Poisson’s ratio
approaches 0.5. On the one hand we show a standard shifted Laplacian approach without using the mixed formulation, applying W(2,2) cycles with damped Jacobi relaxation as preconditioner. On the other hand, we show the shifted Laplacian approach using the mixed formulation (12) using W(1,1) cycles with redblack cellwise relaxation. Table 1 summarizes the experiments. It is clear that the standard shifted Laplacian method performs reasonably well in cases when . Similarly to the acoustic Helmholtz, difficulties are observed when the frequency is high. This behavior is expected. However, the standard method also deteriorates as gets larger ( is fixed), while the proposed approach using the mixed formulation (12) does not deteriorate at all as grows. This is exactly what our method aims to achieve.
3.2 Shifted Laplacian for elastic vs. acoustic Helmholtz
Shifted Laplacian: acoustic (using shear velocity) vs. elastic  

Points per  Acoustic  Elastic  
Grid size  wavelength  #  
15  25  27  
10  40  37  
15  45  47  
10  86  78  
15  75  79  
10  196  148 
Our second set of experiments aims to demonstrate our main message: the new shifted Laplacian method solves the elastic (2) with approximately the same performance of the standard shifted Laplacian for the acoustic equation (1), only with respect to the shear wavenumber instead of the lower pressure wavenumber. To demonstrate this we use the linear velocity and density models in Fig. 4.
We use three gridsizes, for which we solve both the acoustic and elastic equations using 10 and 15 gridpoints per wavelength with respect to the shear velocity. The acoustic equation is defined with the shear velocity instead of the pressure velocity. Both equations are solved with 3level W cycles, and shift parameter . The acoustic equation is solved with cycles with weighted Jacobi as relaxation (weight of 0.8) as in [50]. Table 2 summarizes the results. It is clear that the performance of the shifted Laplacian in the two scenarios is comparable, even though the problems are different, and the discretization and relaxation methods are different. This shows that the main difficulty of the problem lies only in the indefiniteness of the linear systems to be solved.
3.3 Performance comparison of three vs. four levels
Shifted Laplacian: three vs. four levels  
3 levels  3 levels  4 levels  
Grid size  
57  86  93  
84  133  157  
139  193  236 
In this set of experiments we demonstrate the performance of our solver on the highly heterogenous Marmousi2 elastic twodimensional model [33], which appears in Fig. 5. This is a 2D model, but as said before, we consider it as a case study for real 3D scenarios. Since the model is shallow (3 km deep), we extend it by half a km at the bottom to accommodate the absorbing boundary layer. In these experiments we try to address one limitation of our method  the need to exactly solve the coarsest grid problem, which at three levels is still quite large in three dimensions. That is also a problem when solving the acoustic equation, only now the system is much larger (for the same mesh), especially when using the mixed formulation. To reduce the size of the coarsest grid, we inevitably use four levels in our multigrid hierarchy. A similar effort was performed in [10] for the acoustic equation, where the authors show that four levels can be used, but with a lower damping parameter for the relaxation on the third grid. We apply a similar strategy here, using a damping parameter of 0.5 for the Vanka relaxation on the first two grids, and 0.2 on the third. In addition, we also increase the shift parameter to , simply because the 4level approach failed to converge for . The rest of the parameters are the same as the default parameters mentioned earlier.
Table 3 shows the results for the Marmousi2 model, where we choose the frequency to correspond to 10 gridpoints per shear wavelength. In terms of iteration counts it is clear that using 3 levels is much more effective than 4 levels, which is also the conclusion of [10] for the acoustic case. It is also clear that the choice of is significant as the preconditioner using is more effective than the one with . Given , the multigrid cycles with 3 and 4 levels yield comparable performance.
3.4 Threedimensional experiments
In this set of experiments we demonstrate our ability to solve the problem in three dimensions. We use a threedimensional version of the linear model in Fig. 4, and in addition use the heterogenous Overthrust model in Fig. 6. Because the model is shallow, we add 16 grid points at the bottom of the domain to accommodate the absorbing boundary layer. Since this model is only acoustic and includes only pressure velocity, we set the shear velocity to be and set the density to be constant .
The linear systems (12) are huge in 3D, and hence for these experiments we use much smaller grid sizes than in the 2D experiments. We use the fourlevel hierarchy and a shift of , and also use single precision in the multigrid preconditioner computations to save memory. The tests were computed on a workstation with Intel Xeon Gold 5117 2GHz X 2 (14 cores per socket) with 256 GB RAM, running on Centos 7 Linux distribution.
Table 4 summarizes the results. Generally, the iteration counts are comparable to those obtained in the 2D results in Table 3 using four levels. Here again we note that the four levels and the large shift significantly hurt the performance of the solver.
The 3D linear model  The 3D Overthrust model  

Grid size  Grid size  
36  35  
36  50  
79  67  
134  97 
4 Conclusions and future work
In this paper we present a new shifted Laplacian multigrid approach for solving the elastic Helmholtz equation, which has less solvers available in the literature than its acoustic counterpart. The main idea of our approach is that in the case of zero frequency the elastic Helmholtz problem is essentially the linear elasticity problem, which is known to cause difficulties to standard multigrid methods. Hence, our method combines the shifted Laplacian approach with approaches for linear elasticity. That is, with some specialized components in the multigrid cycles, such as mixed formulation and cellwise relaxation, the elastic Helmholtz equation can be solved with the same efficiency as the acoustic Helmholtz equation, only with respect to the shear wavenumber instead of the pressure wavenumber. We verify this property in our experiments, and in addition show that our solver is robust with respect to the Poisson’s ratio  a property that is important in real life scenarios.
In our research we aim to solve the elastic Helmholtz problem in large scale 3D scenarios. However, the linear systems that result from the discretization are huge. Even after coarsening twice, the linear systems (on the third grid) are too large for direct solvers, especially in mixed formulation which does not have a favorable sparsity pattern. In this work we inevitably used four levels in 3D to solve the problem in reasonable resolution using our resources. We showed that this choice is inferior to using three levels.
We aim to solve the 3D problem at larger, more realistic scales, and we need more ways to deal with the problem size. In our future research we will explore ways to efficiently distribute and solve the problem using many computing nodes. In addition, we will explore options to combine the multigrid approach with adaptive discretization, as the shear velocity tend to be much higher in the upper regions of the domain (requiring a finer mesh) than in the lower ones.
Overall, this work demonstrates that besides the size, there is no real “added difficulty” to the elastic Helmholtz equation compared to the acoustic one, and that the two problems have a lot in common. The real question still remains: how to better deal with the indefiniteness of the Helmholtz problems (acoustic or elastic). This is also part of our future research. We think that most approaches for the acoustic equation can be extended to the elastic equation as well.
5 Acknowledgement
The research leading to these results has received funding from the European Union’s  Seventh Framework Programme (FP7/20072013) under grant agreement number 623212  MC Multiscale Inversion.
References
 [1] T. Airaksinen, E. Heikkola, A. Pennanen, and J. Toivanen, An algebraic multigrid based shiftedlaplacian preconditioner for the Helmholtz equation, J. Comput. Phys., 226 (2007), pp. 1196–1210.
 [2] P. R. Amestoy, I. S. Duff, J.Y. L’Excellent, and J. Koster, A fully asynchronous multifrontal solver using distributed dynamic scheduling, SIAM Journal on Matrix Analysis and Applications, 23 (2001), pp. 15–41.
 [3] A. Bayliss, C. I. Goldstein, and E. Turkel, On accuracy conditions for the numerical computation of waves, J. Comput. Phys., 59 (1985), pp. 396–404.
 [4] J. Bezanson, A. Edelman, S. Karpinski, and V. B. Shah, Julia: A fresh approach to numerical computing, SIAM Review, 59 (2017), pp. 65–98.
 [5] D. Borisov and S. C. Singh, Threedimensional elastic full waveform inversion in a marine environment using multicomponent oceanbottom cables: a synthetic study, Geophysical Journal International, 201 (2015), pp. 1215–1234.
 [6] A. Brandt and I. Livshits, Waveray multigrid method for standing wave equations, Electron. Trans. Numer. Anal, 6 (1997), pp. 162–181.
 [7] W. L. Briggs, V. E. Henson, and S. F. McCormick, A multigrid tutorial, SIAM, second ed., 2000.
 [8] R. Brossier, Twodimensional frequencydomain viscoelastic full waveform inversion: Parallel algorithms, optimization and performance, Computers & Geosciences, 37 (2011), pp. 444–455.
 [9] R. Brossier, S. Operto, and J. Virieux, Seismic imaging of complex onshore structures by 2d elastic frequencydomain fullwaveform inversion, Geophysics, 74 (2009), pp. WCC105–WCC118.
 [10] H. Calandra, S. Gratton, X. Pinel, and X. Vasseur, An improved twogrid preconditioner for the solution of threedimensional Helmholtz problems in heterogeneous media, Numerical Linear Algebra with Applications, 20 (2013), pp. 663–688.
 [11] W. Chen, Y. Liu, and X. Xu, A robust domain decomposition method for the helmholtz equation with high wave number, ESAIM: Mathematical Modelling and Numerical Analysis, 50 (2016), pp. 921–944.
 [12] R. Clayton and B. Engquist, Absorbing boundary conditions for acoustic and elastic wave equations, Bulletin of the seismological society of America, 67 (1977), pp. 1529–1540.
 [13] S. Cools, P. Ghysels, W. van Aarle, J. Sijbers, and W. Vanroose, A multilevel preconditioned krylov method for the efficient solution of algebraic tomographic reconstruction problems, Journal of Computational and Applied Mathematics, 283 (2015), pp. 1–16.
 [14] S. Cools, B. Reps, and W. Vanroose, A new leveldependent coarse grid correction scheme for indefinite Helmholtz problems, Numerical Linear Algebra with Applications, 21 (2014), pp. 513–533.
 [15] H. C. Elman, O. G. Ernst, and D. P. O’leary, A multigrid method enhanced by krylov subspace iteration for discrete helmholtz equations, SIAM J. Sci. Comput., 23 (2001), pp. 1291–1315.
 [16] B. Engquist and A. Majda, Radiation boundary conditions for acoustic and elastic wave calculations, Communications on pure and applied mathematics, 32 (1979), pp. 313–357.
 [17] B. Engquist and L. Ying, Sweeping preconditioner for the Helmholtz equation: hierarchical matrix representation, Communications on pure and applied mathematics, 64 (2011), pp. 697–735.
 [18] I. Epanomeritakis, V. Akcelik, O. Ghattas, and J. Bielak, A NewtonCG method for largescale threedimensional elastic fullwaveform seismic inversion, Inverse Problems, (2008).
 [19] Y. A. Erlangga and R. Nabben, On a multilevel Krylov method for the Helmholtz equation preconditioned by shifted laplacian, Electronic Transactions on Numerical Analysis, 31 (2008), p. 3.
 [20] Y. A. Erlangga, C. W. Oosterlee, and C. Vuik, A novel multigrid based preconditioner for heterogeneous Helmholtz problems, SIAM J. Sci. Comput., 27 (2006), pp. 1471–1492.
 [21] M. J. Gander and H. Zhang, Domain decomposition methods for the helmholtz equation: a numerical investigation, in Domain Decomposition Methods in Science and Engineering XX, Springer, 2013, pp. 215–222.
 [22] F. Gaspar, J. Gracia, F. Lisbona, and C. Oosterlee, Distributive smoothers in multigrid for problems with dominating grad–div operators, Numerical linear algebra with applications, 15 (2008), pp. 661–683.
 [23] D. Gordon and R. Gordon, Robust and highly scalable parallel solution of the Helmholtz equation with large wave numbers, Journal of Computational and Applied Mathematics, 237 (2013), pp. 182–196.
 [24] B. GosselinCliche and B. Giroux, 3d frequencydomain finitedifference viscoelasticwave modeling using weighted average 27point operators with optimal coefficients, Geophysics, 79 (2014), pp. T169–T188.
 [25] E. Haber and S. MacLachlan, A fast method for the solution of the Helmholtz equation, J. Comput. Phys., 230 (2011), pp. 4403–4418.
 [26] K. Kelly, R. Ward, S. Treitel, and R. Alford, Synthetic seismograms: A finitedifference approach, Geophysics, 41 (1976), pp. 2–27.
 [27] J. Kormann, J. E. Rodríguez, M. Ferrer, A. Farrés, N. Gutiérrez, J. de la Puente, M. Hanzich, and J. M. Cela, Acceleration strategies for elastic full waveform inversion workflows in 2d and 3d, Computational Geosciences, 21 (2017), pp. 31–45.
 [28] A. R. Levander, Fourthorder finitedifference psv seismograms, Geophysics, 53 (1988), pp. 1425–1436.
 [29] Y. Li, B. Han, L. Métivier, and R. Brossier, Optimal fourthorder staggeredgrid finitedifference scheme for 3D frequencydomain viscoelastic wave modeling, J. Comput. Phys., 321 (2016), pp. 1055–1078.
 [30] Y. Li, L. Métivier, R. Brossier, B. Han, and J. Virieux, 2D and 3D frequencydomain elastic wave modeling in complex media with a parallel iterative solver, Geophysics, 80 (2015), pp. T101–T118.
 [31] Q. Liao and G. A. McMechan, Multifrequency viscoacoustic modeling and inversion, Geophysics, 61 (1996), pp. 1371–1378.
 [32] I. Livshits, A scalable multigrid method for solving indefinite Helmholtz equations with constant wave numbers, Numerical Linear Algebra with Applications, 21 (2014), pp. 177–193.
 [33] G. S. Martin, R. Wiley, and K. J. Marfurt, Marmousi2: An elastic upgrade for marmousi, The Leading Edge, 25 (2006), pp. 156–166.
 [34] L. Métivier, R. Brossier, S. Operto, and J. Virieux, Full waveform inversion and the truncated newton method, SIAM Review, 59 (2017), pp. 153–195.
 [35] L. N. Olson and J. B. Schroder, Smoothed aggregation for Helmholtz problems, Numerical Linear Algebra with Applications, 17 (2010), pp. 361–386.
 [36] C. Oosterlee, C. Vuik, W. Mulder, and R.E. Plessix, Shiftedlaplacian preconditioners for heterogeneous Helmholtz problems, in Advanced Computational Methods in Science and Engineering, Springer, 2010, pp. 21–46.
 [37] S. Operto, J. Virieux, P. Amestoy, J.Y. L’Excellent, L. Giraud, and H. B. H. Ali, 3d finitedifference frequencydomain modeling of viscoacoustic wave propagation using a massively parallel direct solver: A feasibility study, Geophysics, 72 (2007), pp. SM195–SM211.
 [38] J. Poulson, B. Engquist, S. Li, and L. Ying, A parallel sweeping preconditioner for heterogeneous 3D Helmholtz equations, SIAM J. Sci. Comput., 35 (2013), pp. C194–C212.
 [39] R. Pratt, Seismic waveform inversion in the frequency domain, part 1: Theory, and verification in a physical scale model, Geophysics, 64 (1999), pp. 888–901.
 [40] B. Reps and T. Weinzierl, Complex additive geometric multilevel solvers for helmholtz equations on spacetrees, ACM Transactions on Mathematical Software (TOMS), 44 (2017), p. 2.
 [41] G. Rizzuti and W. A. Mulder, Multigridbased shiftedlaplacian preconditioning for the timeharmonic elastic wave equation, J. Comput. Phys., 317 (2016), pp. 47–65.
 [42] L. Ruthotto, E. Treister, and E. Haber, jInv – a flexible Julia package for PDE parameter estimation, SIAM J. Sci. Comput., 39 (2017), pp. S702––S722.
 [43] Y. Saad, A flexible innerouter preconditioned GMRES algorithm, SIAM J. Sci. Comput., 14 (1993), pp. 461–469.
 [44] O. Schenk and K. Gärtner, Solving unsymmetric sparse systems of linear equations with pardiso, Future Generation Computer Systems, 20 (2004), pp. 475–487.
 [45] I. Singer and E. Turkel, Highorder finite difference methods for the helmholtz equation, Computer Methods in Applied Mechanics and Engineering, 163 (1998), pp. 343–358.
 [46] I. Singer and E. Turkel, A perfectly matched layer for the helmholtz equation in a semiinfinite strip, J. Comput. Phys., 201 (2004), pp. 439–465.
 [47] I. Singer and E. Turkel, Sixthorder accurate finite difference schemes for the helmholtz equation, Journal of Computational Acoustics, 14 (2006), pp. 339–351.
 [48] I. Štekl and R. G. Pratt, Accurate viscoelastic modeling by frequencydomain finite differences using rotated operators, Geophysics, 63 (1998), pp. 1779–1794.
 [49] C. C. Stolk, A rapidly converging domain decomposition method for the helmholtz equation, J. Comput. Phys., 241 (2013), pp. 240–252.
 [50] E. Treister and E. Haber, Full waveform inversion guided by travel time tomography, SIAM J. Sci. Comput., 39 (2017), pp. S587––S609.
 [51] , A multigrid solution to the helmholtz equation with a point source based on traveltime and amplitude, Submitted, (2018).
 [52] U. Trottenberg, C. Oosterlee, and A. Schüller, Multigrid, Academic Press, London and San Diego, 2001.
 [53] P. Tsuji and R. Tuminaro, Augmented amgshifted laplacian preconditioners for indefinite helmholtz problems, Numerical Linear Algebra with Applications, 22 (2015), pp. 1077–1101.
 [54] E. Turkel, D. Gordon, R. Gordon, and S. Tsynkov, Compact 2d and 3d sixth order schemes for the helmholtz equation with variable wave number, J. Comput. Phys., 232 (2013), pp. 272–287.
 [55] N. Umetani, S. P. MacLachlan, and C. W. Oosterlee, A multigridbased shifted laplacian preconditioner for a fourthorder helmholtz discretization, Numerical Linear Algebra with Applications, 16 (2009), pp. 603–626.
 [56] H. A. Van der Vorst, BiCGSTAB: A fast and smoothly converging variant of BiCG for the solution of nonsymmetric linear systems, SIAM Journal on scientific and Statistical Computing, 13 (1992), pp. 631–644.
 [57] P. Vanek, J. Mandel, and M. Brezina, Algebraic multigrid by smoothed aggregation for second and fourth order elliptic problems, Computing, 56 (1996), pp. 179–196.
 [58] S. Vanka, Blockimplicit multigrid calculation of twodimensional recirculating flows, Computer Methods in Applied Mechanics and Engineering, 59 (1986), pp. 29–48.
 [59] J. Virieux, Psv wave propagation in heterogeneous media: Velocitystress finitedifference method, Geophysics, 51 (1986), pp. 889–901.
 [60] J. Virieux and S. Operto, An overview of fullwaveform inversion in exploration geophysics, Geophysics, 74 (2009), pp. WCC1–WCC26.
 [61] H. Wobker and S. Turek, Numerical studies of vankatype smoothers in computational solid mechanics, Advances in Applied Mathematics and Mechanics, 1 (2009), pp. 29–55.
 [62] Y. Zhu, E. Sifakis, J. Teran, and A. Brandt, An efficient multigrid method for the simulation of highresolution elastic solids, ACM Transactions on Graphics (TOG), 29 (2010), p. 16.