A damping preconditioner for time-harmonic wave equations in ﬂuid and elastic material

A physical damping is considered as a preconditioning technique for acoustic and elastic wave scattering. The earlier preconditioners for the Helmholtz equation are generalized for elastic materials and three-dimensional domains. An algebraic multigrid method is used in approximating the inverse of damped operators. Several numerical experiments demonstrate the behavior of the method in complicated two-dimensional and three-dimensional domains


Introduction
Developing efficient methods to solve acoustic and elastic scattering problems has proved to be challenging by mathematical and computational means.These problems have a wide range of applications in different disciplines, and therefore there is a big interest to find efficient methods to solve these problems numerically.Modeling is done by acoustic or elastic wave equation, depending on the material, and it is often sufficient to consider only time-harmonic solutions.For incompressible fluids, the reduced wave equation is the Helmholtz equation.For linearly elastic material, the Navier equation can be applied.An approximate solution can be obtained by discretizing these equations using, for example, a finite difference or finite element method.
Finite element methods have become a popular technique to discretize partial differential equations in complex geometries.It has successfully been used for interior scattering problems like acoustic scattering in a car cabin [1] as well as for exterior problems.A review [2] gives an overview of recent research on finite element methods for acoustic problems.Since the paper [3] the research on the construction of absorbing boundary conditions and absorbing layers at the truncation boundary of the exterior domain has been active; see [2] and references therein.The size of the scattering problems is often limited in high frequency problems because the methods become ineffective as the frequency grows.Particularly the finite element phase shift (pollution) error necessitates finer discretizations for high-frequency problems [4] and thus an increasing memory and computational requirements.
The resulting systems of linear equations from the discretization of the Helmholtz equation and the Navier equation are non-Hermitian and indefinite, and for mid-frequency and high-frequency problems, they can be extremely large.These properties make them a challenge for the current solvers.For two-dimensional problems, it is often feasible to use direct methods for solving these systems, but three-dimensional problems lead to systems that can not be solved by these methods with an affordable computing effort.Hence, it is necessary to use iterative methods such as the GMRES method [5] or the Bi-CGSTAB method [6].However, these methods require a good preconditioner for the discretized equations in order to have reasonably fast convergence.
Several preconditioners and iterative solution techniques have been proposed for the discrete Helmholtz and Navier equations.Domain decomposition methods have been proposed for Helmholtz problems in [7,8,9,10,11,12], and for elastic problems in [13,14,15,16].Controllability methods have been proposed for both Helmholtz and Navier problems in [17,18].Multigrid methods have been considered for acoustic and elastic problems in [19,20,21,22].With multigrid methods, it is difficult to define a stable and sufficiently accurate coarse grid problem and smoother for it.For acoustic and elastic problems in homogenous medium, domain imbedding/fictitious domain methods in [23,24,25,26] have been fairly effective, but these methods are pretty restrictive and not well-suited in general, complicated domains.An incomplete factorization preconditioner has been considered in [27], for example, and in [28] a tensor product preconditioner is used.
So called natural preconditioning techniques are applicable for many problems including time-harmonic wave equations [29].The class of preconditioners based on damped operators that are considered here, is an example of this approach.A shifted-Laplacian preconditioner with a complex shift, which is called here a damped Helmholtz preconditioner, was first considered in [30] for the Helmholtz equation.This was a development over the shifted-Laplacian preconditioner with a real shift previously described in [31].Already in [32,33] a complex shift was employed, but for a completely different way and purpose: it was used to transform a singular problem into a non singular one.Here the purpose to introduce a complex shift into a preconditioner for a non singular problem is to enable the effecient use of multigrid methods.
A damped Helmholtz preconditioner with geometric multigrid was considered in [21].There, the scattering problems were posed in a rectangular domain and they were discretized using low-order finite differences.Our earlier study [34] extended this approach to general shaped two dimensional domains using linear, quadratic, and cubic finite element discretizations by applying an algebraic multigrid (AMG) instead of the geometric multigrid to approximate the inversion of the damped Helmholtz operator.In [35], this method was compared with the previously mentioned controllability method.
In this paper, a generalization will be proposed to the preconditioner described in [34], an AMG-based damped preconditioner for time-harmonic wave propagation problems in elastic media, i.e. the Navier equation.This preconditioner will be called a damped Navier preconditioner.Results considering the eigenvalue spectrum of the shifted-Laplacian preconditioned discretized Helmholtz equation were given in [36] and some of these will be generalized to the Navier equation.Simulations are carried out in twodimensional and three-dimensional computational domains including complicated geometries for both Helmholtz and Navier problems.This paper is organized as follows.In Section 2 acoustic and elastic wave scattering models and their discretizations are described.The iterative solution and preconditioning are discussed in Section 3 and mathematical results on the eigenvalue spectrum are given in Section 4. The algebraic multigrid method employed in the preconditioning is described in Section 5. Then numerical results are presented in Section 6 and finally, conclusions are given in Section 7.

Wave scattering in fluids
For a time-harmonic pressure of the form p (x, t) = p (x) e −iωt with an an- gular frequency ω and imaginary unit i = √ −1, the wave scattering in a fluid domain Ω f can be described by a Helmholtz equation where and where g f (x) describes a sound source and n(x) is the outer normal vector.Choosing the absorbency coefficient γ to be zero leads to the Neumann boundary condition and γ = 1 gives a low-order absorbing boundary condition.

Wave scattering in elastic materials
For time-harmonic displacements u (x, t) = e −iωt û (x) in a domain Ω s con- sisting of elastic materials, the scattering of time-harmonic waves can be described by a Navier equation where σ is the stress tensor, f s is a force term, and ρ s (x) is the density of the material.Hooke's law gives a relation between displacements, and stress and strain forces, thus describing strain tensor ǫ and stress tensor σ: Here Lamé parameters λ and µ are defined as follows: These depend on Young modulus E (x) and the Poisson ratio ν (x) that char- acterize the elastic behaviour of the material.The speed of pressure wave, c p , and shear wave, c s , can be expressed as functions of Lamé parameters: Wavelengths and wave numbers for pressure and shear waves are where g s (x) describes the vibration source.The impedance boundary con- dition on Γ s i is approximated by the equation where γ is the absorbency coefficient, i = √ −1, B is a 2 × 2 matrix for twodimensional problems (D = 2) and a 3 × 3 matrix for three-dimensional problems (D = 3).Choosing the absorbency coefficient γ = 0 leads to natural boundary condition and γ = 1 gives an absorbing boundary condition.In component form B has expressions for D = 2 and where c p and c s are the speeds of pressure and shear waves given by (7) and n = (n 1 , • • • , n D ) T is the normal vector pointing out of elastic domain, and t = (t 1 , • • • , t D ) T and s = (s 1 , • • • , s D ) T are tangential vectors on the boundary.

Weak formulation and finite element discretization
For the weak formulation of the Helmholtz equation, we define a test function space V f 0 and a solution space V f g as The weak form of (1) reads: Find p ∈ V f g such that for all q ∈ V f 0 .Similarly, for the Navier equation, we define a test function space V s 0 and a solution space V s g as Now, the weak form of (4) reads: Find û ∈ V s g such that for all v ∈ V s 0 .
For a finite element discretization, a mesh K h is defined such that Ω h = τ∈K h τ.The mesh consists of triangles τ in two-dimensional and of tetrahedra in three-dimensional problems.Here h denotes the diameter of the largest triangle or tetrahedron and Ω f ,s h is an approximation of Ω f ,s .For the finite elements of order m discrete test function spaces are 16) where P m denotes polynomials of order m.Discrete solution spaces V f ,s g,h are the same except the zero boundary value on Γ f ,s d,h is replaced by approximations of g f and g s .In this paper, linear, quadratic, and cubic finite elements are employed, i.e. m = 1, 2, or 3.For the spaces V h and V g,h , Lagrangian polynomials are used as basis functions.
For the analytical study of eigenvalue spectra in Section 4, it is practical to define the following matrices based on the integrals in ( 13) and (15): [36], the discretized Helmholtz and Navier operators have matrix forms respectively.Now, let the vector w contain the nodal values of p or û, so that for the Helmholtz problem it has form w = [ p1 , • • • , pn ] T , and for the two-dimensional Navier problem it has form . By replacing the spaces, domains, and boundaries in ( 13) or ( 15) by their discrete counterparts, the system of linear equations is obtained.The complex-value sparse matrix A is given by F or S in (18) and f is a vector resulting from an inhomogeneous Dirichlet boundary value and/or a non-zero f f in (13) or f s in (15).
The approximation properties of such finite element discretizations for the Helmholtz equation have been studied in [4].Due to the pollution (phase shift) error, a non-optimal L 2 error estimate is obtained, where C 1 and C 2 are constants.Based on this estimate, larger mesh step sizes can be used when higher order finite elements are being used, in order to reach the same accuracy level.

Iterative solution and damped preconditioner
The matrix A in ( 20) is indefinite and symmetric, but not Hermitian.For example, the generalized minimal residual (GMRES) method [5] and the Bi-CGSTAB method [6] are suitable iterative methods for these equations.
These and other applicable iterative methods are described in [37].The GM-RES method minimizes the 2-norm of the residual on Krylov subspaces.This is a desirable property leading to a monotonic reduction of the norm of the residual over iterations, but a disadvantage is that all basis vectors for the Krylov subspace needs to be stored.This makes the computational cost of the GMRES methods grow quadratically with iterations and also causes linear growth in memory requirement.The computational cost of the Bi-CGSTAB method grows linearly with the iterations and the memory requirement is constant, but the convergence can be erratic and slower than with the GMRES method.In the numerical experiments, the full GMRES method is used without restarts.
The convergence of Krylov subspace methods for the system ( 20) is very slow for medium-and large-scale scattering problems due to the ill conditioning of A. To improve the conditioning and the speed of convergence, a right preconditioner denoted by B is introduced.This leads to a preconditioned system Once ũ is solved from this system, the solution u is obtained as u = B −1 ũ.The goal is to find such a preconditioner B that the matrix AB −1 is well conditioned and that vectors can be multiplied by B −1 , i.e. solve systems with B with a small computational effort.These properties would lead to a fast convergence of the iterative method and to a small overall computational cost.
A shifted-Laplacian with a complex shift z 2 = α 2 + β 2 i was suggested in [30] as a preconditioner for the Helmholtz equation.By choosing α 2 = 1 and β 2 to be negative, F d is the Helmholtz operator in (1) with some additional damping.Using the matrices defined in (18), the discretization of F d leads to a matrix With sufficient damping, systems with F d can be solved much more easily than with F and the conditioning of FF −1 d can still be good.The use of different approximations for F −1  d have been studied in [30,21,38,34].Here an algebraic multigrid approximation described in Section 5 is considered.
Our hypotesis is that a similar physical damping can be employed to construct an efficient preconditioner for the Navier equations.Damping in elastic materials can be modelled by using a complex Young modulus.Multiplying the original Young modulus E(x) by a complex z 2 leads to a precon-ditioning operator The coefficient z 2 appears also in the impedance boundary condition (10) as follows Using the matrices in (18), the discretization of S d leads to the damped Navier preconditioner 4 Spectral analysis for the preconditioned Navier equation Studying the eigenvalue spectrum of the preconditioned matrix AB −1 is an usual way to estimate the convergence of an iterative method like GMRES.
In [36], Theorems 3.1-3.6give useful information of the eigenvalue spectrum of the preconditioned Helmholtz operator.Some of these results can be generalized to the Navier equation, as will be shown in the following.
As defined in (19), the matrix of the discretized Navier equation is Here matrices L s and C s are symmetric positive semi-definite and M s is symmetric positive definite, and z 1 is a complex number.The case that there are only natural and/or Dirichlet boundary conditions, i.e.C s = 0, and the material is not absorbing, is analyzed first.Thus, the matrix S simplifies to The eigenvalue problem AB −1 ỹ = τ ỹ is equivalent to where y = B −1 ỹ.From this, the eigenvalue problem can be derived.
As the matrix L s is positive semi-definite and M s is symmetric positive definite, the eigenvalues λ are real.The eigenvalue τ is a function of λ given by By the change of variable λ ′ = λ −1 , the form is obtained.This is the same equation of a circle in the complex plane that was found in [36] for the eigenvalue spectrum of the preconditioned Helmholtz equation.Due to this, the following corollary of Theorems 3.1-3.3 in [36] can be formulated.
Corollary 1 For the eigenvalues τ = τ r + iτ i of the generalized eigenvalue problem Sy = τS d y, the following statements hold: • If β 2 = 0, the eigenvalues are located on straight line in the complex plane given by the equation • If β 2 = 0, the eigenvalues are located in complex plane on the circle given by The center of the circle is at c = z 1 −z 2 z 2 −z 2 and the radius is the origin is not enclosed by the circle defined by (35).
The case of impedance boundary conditions with γ = 0 in (10), i.e.C s = 0, is considered next.It is not evident that the results presented in [36] for the Helmholtz equation with C f = 0 are applicable for the Navier equation.However, numerical experiments in Section 6 suggest that similar behaviour to the one described by Theorems 3.4-3.6 in [36] holds also here.The following states this as a conjecture.
Conjecture 2 For the eigenvalues τ = τ r + iτ i of the generalized eigenvalue problem Sy = τS d y, the following statements hold: • If β 2 = 0, the eigenvalues are located in the half-plane • If β 2 > 0 the eigenvalues are inside or on the circle with the center at c = z 1 −z 2 z 2 −z 2 and the radius R = z 2 −z 1 z 2 −z 2 .If β 2 < 0, the eigenvalues are outside or on the same circle.

Algebraic multigrid based damped preconditioners
The approximation of the inverse of the damped operator B −1 given by a multigrid method is denoted by B −1 MG .In [21], Erlangga, Oosterlee, and Vuik used one cycle of a geometric multigrid method for this.For low-frequency problems the conditioning of AB −1  MG is good.For high-frequency problems the conditioning deteriorates so that the number of Bi-CGSTAB iterations appeares to grow linearly with frequency in [21].They also showed that this preconditioner is well-suited for problems with a varying speed of sound.In [34], the geometric multigrid method was replaced by a more generic and more flexible algebraic multigrid method (AMG).In this paper, an AMG based on [39] is utilized, using the implementation that is described in [34], with modifications that make it suitable for vector valued problems, like the Navier equation.
The employed AMG method uses a graph to construct coarse spaces.Here the graph is based on the discretization mesh.Alternative approach would be to build the graph based on the matrix B. When using linear elements in a scalar problem, both approaches result in the same graph.For an elastic solid modelled by the Navier equation, the graph is formed without connections (edges) between displacement components.This choice is made for two reasons: Adding these connections would cause too rapid coarsening process.Secondly, the error behaves smoothly for each component separately and the AMG method is especially efficient at reducing smooth error components.The graph therefore consists of separate disconnected graphs, one for each displacement component.
For linear finite elements, the initial graph G 0 is the graph defined by the triangulation.For quadratic and cubic elements, the graph is defined by a refined mesh.In two-dimensional domains, quadratic triangle elements are divided into four triangles by connecting the midpoints of the edges, and cubic triangle elements are divided into nine triangles.In three-dimensional domains, quadratic tetrahedron elements are divided into eight and cubic tetrahedra into 26 tetrahedra.If the graph defined by B was used directly with high-order elements, the coarsening procedure would coarsen the graph too rapidly, leading to an impaired conditioning of AB −1  MG and a slower convergence of the GMRES method.
The nodes onto a coarser graph G k+1 are chosen from the nodes of G k as follows.Find the node in G k which has the smallest degree, i.e. the smallest number of edges associated to it.If there are several such nodes, choose the first one according to the node numbering.This node is included onto the graph G k+1 .Eliminate this node and all its neighbors from the graph G k .Repeat this procedure until there are no nodes left in G k .After choosing the nodes on G k+1 , they are numbered following their order in the numbering of the nodes on G k .
On coarse levels, different displacement components are chosen to be disconnected.Thus, the restriction matrix is defined blockwise as The elements of the diagonal blocks of the restriction matrix are defined by the rule 1 for a fine node j which is a coarse node i, for a fine node j which is a neighbor of coarse node i and has n neighboring coarse nodes, 0 otherwise, where fine and coarse refers to the graphs G k and G k+1 , respectively.The edges of the coarse graph G k+1 are formed using the restriction matrix R k .Each coarse graph node corresponds to a row in the restriction matrix.There is an edge between two nodes if and only if the corresponding rows of the restriction matrix have a non-zero element in the same column.
The coarse level matrices are now defined as follows The usual multigrid W-cycle is used with the AMG method.For preconditioning, the initial approximate solution is zero in the multigrid algorithm.At each level, presmoothing and postsmoothing is performed by one underrelaxed Jacobi iteration.At the coarsest level, a direct solver is used instead of an iterative method.

Numerical results
Numerical simulations were carried out on selected example problems.In Subsection 6.1, the eigenvalues of two-dimensional Navier problems are studied and compared with the results presented in Section 4. In Subsection 6.2, the performance of the method is considered for two-dimensional and three-dimensional Helmholtz and Navier problems by measuring iteration counts required to satisfy a convergence criterion.
The following material parameters are used in tests unless specified otherwise.The Helmholtz problems have domain Ω f consisting of air, with the density ρ f = 1.2 kg/m 3 and the speed of sound c = 344 m/s.The Navier problems are posed in a domain Ω s consisting of aluminum with the density ρ s = 2700 kg/m 3 , Young modulus E = 7.00 • 10 10 Pa, and Poisson ratio ν = 0.33.Meshes were generated using Comsol Multiphysics 3.3 in such a way that the maximum element size is h = λ/10, where λ is the wavelength of slowest wave mode.In the Helmholtz problems, λ is the length of acoustic waves, and in the Navier problems, it is the length of shear waves.

Eigenvalues
In [34], the eigenvalue spectra of the preconditioned system matrices were examined for several two-dimensional Helmholtz example problems.Here the eigenvalue spectra will be studied, when the system is preconditioned by a damped preconditioner for two-dimensional and three-dimensional Helmholtz and Navier problems.Two-dimensional problems are studied in the unit square domain like in [21,34] for the Helmholtz problem.A three-dimensional cube domain will also be considered for both Helmholtz and Navier problems.Estimates for the eigenvalue spectra of the preconditioned Navier equation, when Dirichlet or absorbing boundary conditions are posed on boundaries were presented in Section 4. These estimates will be compared to the numerically obtained eigenvalues.
First, the unit square problem will be considered for the Navier equation.The frequency 2.2 kHz is used in the eigenvalue study.The eigenvalues of AB −1 for the unit square problem with the Dirichlet and absorbing boundary conditions are presented in Figure 1.Also the eigenvalues of AB −1 MG are plotted for the same problem, where B −1 MG is the algebraic multigrid approximation of B −1 .The eigenvalue spectrum for the Navier problem with Dirichlet boundary conditions is distributed exactly on the circle as (35) describes.It is also seen that the algebraic multigrid does not spread the spectrum much.Most of eigenvalues seem to move slightly closer to the center of the circle.For the eigenvalue spectra of the problems with the absorbing boundary conditions, it will be shown that the inequality (36) holds in numerical examples.Similar inequality was proven in [36] to hold for Helmholtz problems.According to the inequality (36), the eigenvalues should lie inside or outside of the circle depending on the sign of β 2 .For better conditioning, β 2 is always chosen positive.Thus, according to (36), the eigenvalues should lie inside the circle.For the unit square problem with the absorbing boundary conditions, the conjecture seems to be valid, as can be seen in Figure 1.The algebraic multigrid approximation changes the spectrum, but the eigenvalues seem to still lie inside the circle.
For three-dimensional experiments, the cube (0.3 m) 3 is discretized by us- ing linear finite elements for both Helmholtz and Navier problems.For the Navier problem, the frequency f is 5 kHz and for the Helmholtz problem, the frequency f is 500 Hz.In Figure 2, the eigenvalues of AB −1 are plotted and in Figure 3, the eigenvalues of the system with the AMG approximation of the inverse of the damped operator, AB −1  MG , are plotted.Also the circle (35) is drawn in these figures.It is clearly seen in Figure 2, that both Corollary 1 and Conjecture 2 holds for this problem.

Performance of the preconditioner
The performance of the damped preconditioner with the algebraic multigrid will be reported for several different test problems.First, a two-dimensional Navier problem is studied in the unit square and three-dimensional Helmholtz and Navier problems are studied in a cube domain.Then, the method is tested on complicated three-dimensional problems: For the Helmholtz equation, a three-dimensional car cabin domain and a layered wedge domain with a varying speed of sound are considered.For the Navier equation, a crankshaft geometry defined by a Comsol Multiphysics 3.3 example problem is considered.The iteration counts give the number of iterations needed to reduce the relative residual to 10 −6 .
In all performance studies with the unit square problems, the absorbing boundary condition given by (10) with γ = 1 was posed on the boundaries.For the Navier equation, the best value for β 2 was determined as follows.

Cube problem
The Helmholtz and Navier problems were solved in the cube (0.3 m) 3 with a point source in the middle.The performance of the damped preconditioner was compared to a modified incomplete Cholesky factorization (MIC) preconditioner [27].The algorithm presented in [40] is used for the MIC(l) approximation of A −1 , where the parameter l describes the level of fillin in the factorization.The values l = 0 and 1 have been used as bigger ls were uncompetitive as forming incomplete factorization required much more computation.The performance was measured by the number of GM-RES iterations and the total number of floating point operations (FLOPs) required by the preconditioning.The number of FLOPs is a good measure, as it includes the initialization process in addition to GMRES iteration.The results are presented in Tables 2 and 3.
MIC(1) seems to require fewer iterations than AMG, whereas MIC(0) requires more iterations.Especially with linear finite elements, the convergence with MIC(1) is faster than with AMG and MIC(0), as can be seen in Table 2.However, Table 3 shows that the number of FLOPs with the MIC(1) preconditioner is about twice the number with AMG.This is mainly due to expensive factorization process before the iteration.With quadratic elements the MIC preconditioner seems to perform much worse than AMG, both in iteration counts as well as FLOPs.This is true for both Navier and Helmholtz problems.MIC(1) is also using more memory than AMG, although the difference is not substantial.The car cabin problem is a three-dimensional generalization of the twodimensional car cabin problem in [34].The sound source is modelled as the Dirichlet boundary condition p = 1 posed on the wall behind pedals.The impedance boundary condition (3) with γ = 0.2 is posed on the other boundaries.The height of the car cabin is 1.5 meters, the width is 1.5 meters, and the length is 3 meters.An example solution is plotted in Figure 5. Iteration counts are reported in Table 4.For this problem also, the number of iterations grow roughly linearly with respect to the frequency.The three-dimensional wedge problem [41] in the unit cube [0, 1] 3 is a generalization of a two-dimensional problem studied in [28,21,34].In this acous-tic scattering problem, the material is inhomogeneous leading to a piecewise constant speed of sound.The domain has three layers separated by two planes defined by the equations z = 0.1x + 0.2y + 0.6 and z = −0.2x− 0.15y + 0.4.The speeds of sound from the top layer to the bottom layer are c 1 = 1, c 2 = 1 2 , and c 3 = 5 6 .A point source is placed at (0.5, 0, 0.5) and the absorbing boundary conditions are posed on the boundaries.In the AMG method, the Jacobi relaxation parameter is ω = 0.3.

Three-dimensional car cabin problem
The solution of the Helmholtz equation at f = 2.5 Hz is shown in Figure 6.The iteration counts are reported in Table 5.The same linear growth can be observed as in the previous problems.  .The right end is rigid, i.e. the Dirichlet boundary condition u = (0, 0, 0) is posed on there.Other bound- aries have natural boundary conditions, i.e. the impedance boundary condition (10) with γ = 0.The mesh is made of quadratic finite elements and 300, 000 nodes.
The GMRES iteration counts on a range of frequencies are given in Table 6.As there is no absorption, the problem is singular at some frequencies.Due to this the iteration counts do not behave linearly with respect to the frequency.The solution at f = 3 kHz is illustrated in Figure 7.

Conclusions
A damped Navier preconditioner based on an algebraic multigrid method was introduced for time-harmonic elasticity problems.This is a generalization of a shifted-Laplacian preconditioner for the Helmholtz equation.These preconditioners are efficient for Helmholtz and Navier problems in complicated two-dimensional and three-dimensional domains.Higher-order finite elements can be used for the discretization and Helmholtz problems can have variable coefficients.The proposed approach is especially wellsuited for low-frequency and mid-frequency problems.For high frequencies, iteration counts grow roughly linearly with respect to the frequency.The same behavior was also observed in [34,21].
The performance was compared to a modified incomplete Cholesky (MIC) preconditioner.The AMG based damped preconditioner was more efficient as its initialization requires much less computations than the expensive incomplete factorization procedure.Especially with quadratic finite elements, the AMG preconditioner was clearly faster.
The eigenvalues of the preconditioned system were also studied.The earlier results for the Helmholtz equation in [36] were generalized for the Navier equation.It was shown that the eigenvalues of the preconditioned system with the damped Navier preconditioner are on a circle in the complex plane for Dirichlet and Neumann boundary value problems.When one algebraic multigrid cycle is used instead of the exact inverse of the damped Navier operator, the eigenvalues are spread to some extent, but the conditioning is still fairly good.

Figure 1 .
Figure 1.Eigenvalue plots for the Navier problem in the unit square.The upper plots are for Dirichlet boundary value problems and on the lower plots, the absorbing boundary conditions are posed.On the left plots, the eigenvalues of AB −1 are shown.On the right plots, the eigenvalues of AB −1MG are shown.The circle is defined by(35).The damping parameters are z 1 = 1.0 and z 2 = 1.3 + 0.7i.

Figure 2 .
Figure 2. The eigenvalues of AB −1 .The upper plots are for the Helmholtz problems and the lower ones are for the Navier problems.The left plots are for the Dirichlet boundary value problems and the right ones are for problems with absorbing boundary conditions.The damping parameters are z 1 = 1.0 and z 2 = 1.0 + 0.5i for the Helmholtz problems, and they are z 1 = 1.0 and z 2 = 1.0 + 0.8i for the Navier problems.

Figure 5 .
Figure 5.The solution of the Helmholtz equation at the frequency f = 880 Hz in the three-dimensional car cabin.

Figure 6 .
Figure 6.The solution of the Helmholtz equation at f = 2.5 Hz for the three-dimensional wedge problem.

Figure 7 .
Figure 7.The propagation of elastic waves in a crankshaft at f = 3 kHz.The color scale indicates the amplitude of the displacement, with blue corresponding to small displacement and red corresponding to large displacement.The Navier equation is solved in a complicated three-dimensional domain defined by the crankshaft model from Comsol Multiphysics 3.3.The length of the crankshaft is 1.0 meter and it is made of structural steel defined by: the density ρ = 7850 kg/m 3 , Young modulus E = 2 • 10 11 Pa, and Poisson ratio ν = 0.33.A tangential vibration source on the left end is given by the Dirichlet boundary condition u = (1, 0, 1).The right end is rigid, i.e. the Dirichlet boundary condition u = (0, 0, 0) is posed on there.Other bound- aries have natural boundary conditions, i.e. the impedance boundary condition(10) with γ = 0.The mesh is made of quadratic finite elements and 300, 000 nodes.
is a time-harmonic sound source and ρ (x) is fluid density.In inhomogeneous medium, the wave number k varies depending on location as the sound speed c varies.The boundary of the fluid domain Ω f is decomposed into a Dirichlet boundary Γ

Table 1
The results for the unit square elasticity problem.Iteration counts are given for linear and quadratic finite elements.The solution of the unit square elasticity problem at frequencies 4.9 kHz, 9.8 kHz, 19.6 kHz, and 39.2 kHz.The absorbing boundary conditions are posed on the boundaries.

Table 2
The iteration counts for the cube problem for the Helmholtz and Navier equations.Some counts are missing as the computations were too demanding.

Table 3
The number of millions of floating point operations (MFLOPs) for the Helmholtz and Navier equations.Some numbers are missing as the computations were too demanding.

Table 4
The number of iterations for the three-dimensional car cabin problem for the Helmholtz equation.
6.2.4 Three-dimensional wedge problem for the Helmholtz equation

Table 5
The iteration counts for the three-dimensional wedge acoustic scattering problem.
6.2.5 Crankshaft vibration problem

Table 6
The GMRES iteration counts for the crankshaft vibration problem.