Nonlinear eigenvalue and frequency response problems in industrial practice

Abstract

We discuss the numerical solution of large scale nonlinear eigenvalue problems and frequency response problems that arise in the analysis, simulation and optimization of acoustic fields. We report about the cooperation with the company SFE in Berlin. We present the challenges in the current industrial problems and the state-of-the-art of current methods. The difficulties that arise with current off-the-shelf methods are discussed and several industrial examples are presented. It is documented that industrial cooperation is by no means a one-way street of transfer from academia to industry but the challenges arising in industrial practice also lead to new mathematical questions which actually change the mathematical theory and methods.

Journal of Mathematics in Industry (2011) 1:7

DOI 10.1186/2190-5983-1-7

R E S E A R C H Open Access

Received: 21 February 2011 / Accepted: 27 July 2011 / Published online: 27 July 2011 2011 Mehrmann, Schrder; licensee Springer. This is an Open Access article distributed under the terms of the Creative Commons Attribution License

Abstract Background: We discuss the numerical solution of large scale nonlinear eigenvalue problems and frequency response problems that arise in the analysis, simulation and optimization of acoustic elds. We report about the cooperation with the company SFE in Berlin. We present the challenges in the current industrial problems and the state-of-the-art of current methods.

Results: The difculties that arise with current off-the-shelf methods are discussed and several industrial examples are presented.

Conclusions: It is documented that industrial cooperation is by no means a one-way street of transfer from academia to industry but the challenges arising in industrial practice also lead to new mathematical questions which actually change the mathematical theory and methods.

Keywords Nonlinear eigenvalue problem frequency response problem complex

symmetric linear system block-Arnoldi method acoustic eld computation

automotive industry

AMS Subject Classication 65F18 15A18

1 Background

1.1 Introduction

Trafc noise emissions by transport vehicles, such as cars, trains or airplanes are one of the key factors restricting the quality of life in urban areas. Acoustic waves in

V Mehrmann ([letter]) C Schrder

Institut fr Mathematik, Ma 4-5, TU Berlin, Strae des 17. Juni 136, 10623 Berlin, Germany e-mail: [email protected]

C Schrdere-mail: [email protected]

Nonlinear eigenvalue and frequency response problems in industrial practice

Volker Mehrmann Christian Schrder

Page 2 of 18 Mehrmann, Schrder

dynamically moving vehicles arise from many different sources, such as, for example, the noise of engines or the vibrations of the structure due to external excitations like road contact or head wind. The reduction of such noise emissions is therefore an important factor in the design and the economic success of new products.

To minimize potential noise emissions already in an early design phase of new products (avoiding the costly production of prototypes), it is necessary to use simulation, optimization and control techniques based on mathematical models of the vehicle. Performing these tasks requires mathematical models that describe the acoustic eld associated with the structure of the complete vehicle including its interaction with the environment such as, for example, the air in and around the vehicle or the road or rail contact. Furthermore, these models must allow to identify and inuence potential noise sources, and as a vision for the future they must allow to minimize the emissions.

Although much research is carried out in universities and industrial research and development departments, the model based minimization of the noise emissions of a complete car or train (not to mention airplane) including the majority of external and internal excitations is still a vision for the future. To achieve this, a joint effort is needed that includes the identication of sources, the construction of adequate mathematical models that incorporate all possible sources for acoustic waves, the analysis of these models, concerning robustness of their descriptions and their potential for optimization techniques, the development of numerical methods for the simulation and optimization of these models, as well as the implementation of these techniques as production software on modern high performance computers for industrial use.

1.1.1 Modeling

The modeling of acoustic elds inside or outside of dynamically moving vehicles typically uses coupled systems of partial differential equations (PDEs), for example, for the generation of noise by vibrating parts, surface contact, engine noise or head wind. These methods are well established in the engineering community using commercial nite element packages [1, 2]. However, the techniques for the resulting systems of PDEs still have many deciencies, in particular the development of appropriate solvers for the solution of the linear and nonlinear systems and eigenvalue problems that have to be solved after discretization. These turn out to be extremely large and ill-conditioned - that is, extremely sensible with respect to perturbations in the problem-dening data - when a reasonably ne three-dimensional model is used; they may consist of hundreds of millions of equations.

An even greater challenge is to use these models and methods within an optimization loop. There is no chance to use classical off-the-shelf optimization methods for these problems, the problem size is just too large. Instead, one currently applies model reduction techniques [3, 4] to approximate the given ne model that is used for the simulation by a rather crude model that is used in the optimization.

To make such an approach viable within a design environment, where not only the geometry and topology of the vehicle and its material parameters are subject to changes, but also the interaction with the environment is rather complicated, it is necessary to develop integrated mathematical techniques for the modeling, simulation,

Journal of Mathematics in Industry (2011) 1:7 Page 3 of 18

and optimization that make as much use as possible of the properties of the underlying physical model and to transfer this into numerical methods and production codes.

1.1.2 Content of the paper

In this paper we will discuss some of the work and some of the challenges in a long-term cooperation with the company SFE GmbH in Berlin, Germany, which produces software for the simulation and optimization of acoustic elds.

The cooperation involves the development (and implementation on current high-performance computers) of numerical methods for the solution of linear systems and nonlinear eigenvalue problems arising from discretized partial differential equations modeling noise emissions of cars and trains.

The paper is organized as follows.

Section 1.2 briey describes the modeling of acoustic elds inside a car as coupling of uid and structure vibrations. In Section 2.1 we study the direct frequency response problem, that is, the numerical solution of a series of large, sparse, complex symmetric linear systems with a varying parameter. We will address in particular in-core/out-of-core storage methods, preconditioning and parallel execution aspects.

Section 2.2 treats modal reduction to reduce the dimension of the generated models. This requires the computation of all eigenvalues of a large sparse nonlinear eigenvalue problems in a given region of the complex plane.

The described problems and results demonstrate many challenges in the transfer of state-of-the-art numerical methods into the industrial practice, and show that there is mutual benet from such a cooperation between industry and academia.

1.2 Mathematical modeling of interior car acoustics

To model the propagation of acoustic waves inside a car, as illustrated in Figure 1, the three-dimensional lossless wave equation is used, which can be obtained from the continuity equation (conservation of mass)

Fig. 1 FEM model of a car structure and illustration of the propagation of acoustic waves inside a car.

t + (

v) = 0, (1)

Page 4 of 18 Mehrmann, Schrder

together with the Euler equations of uid dynamics

Here v = v(x; y; z; t) is the particle velocity,

= (x; y; z; t) is the particle density,

and p = p(x; y; z; t) denotes the pressure, depending on Cartesian coordinates x, y,

z and time t. Under the simplifying assumptions that there is no temperature change, that the uid is inviscid, that is, no shear forces occur, that the inuence of external forces is restricted to those coming from displacements of the structure at the boundaries, that the uid is adiabatic, that is, there is no heat exchange during compression, that we have an ideal gas, that is, = pc2 , where c is the speed of sound, and nally

that (v )v and vt are small, we can make the Taylor expansions

p = p0 + p(x; y; z; t), with p0 p

= 0 + (x; y; z; t), with 0 .

Then from the Euler equation we obtain 0(vt ) = p or by differentiation

0(vt ) = p, and from the continuity equation we get t + 0v = 0. Alto

gether we obtain the highly simplied system of partial differential equations

1 c2

to which we have to add appropriate initial and boundary conditions. In a closed structure like the interior of the car, we can use the boundary conditions that are obtained from the displacement of the structure to obtain the uid structure interaction, in an open structure like the acoustic waves emitted by the car to the outside, we have to incorporate appropriate far eld conditions or transparent boundary conditions [5].

Let u be the vector of displacements of the structure on the surface. Then v = ut

and thus with the outer normal we get

0 2u

t2 = p.

To obtain a variational formulation, we multiply the equations with a test function w, and use integration by parts by integrating over control volumes V with surface elements S. This gives

V1c2 w2pt2 dV +

V (w)p dV = 0 [integraldisplay]Sw 2ut2 dS,

or equivalently

V10c2 w2pt2 dV +

V10 (w)p dV = [integraldisplay]Sw 2ut2 dS.

One of the most difcult tasks in the modeling is the incorporation of appropriate damping and absorption, because this depends in a rather complicated way on the

vt + (v )v

[parenrightbigg] = p.

(2)

p0 = 106p

,

2p

t2 + 0

v

t = p + 0

v

t = 0

Journal of Mathematics in Industry (2011) 1:7 Page 5 of 18

materials used inside the car, the surface geometry of the interior and this is also one of the terms, where the inuence of the optimization is strong.

In the model treated here, damping and absorption was realized by an additional term

S wr20c2pt dS,

where r = r() is depending on the damping properties of the material (described by

a parameter vector ). We then obtain

V10c2 w2t2 p dV +

S wr20c2pt dS +

V10 (w)p dV = [integraldisplay]Sw 2ut2 dS.

To this variational formulation we can directly apply nite element discretization in space. This leads (in time) to a second order system of implicit differential equations in descriptor form

Mf pd + Df pd + Kf pd + Dsfd = 0, (3)

where Mf = MTf is a positive denite mass matrix, Kf = KTf is a positive denite

stiffness matrix, Df () is a symmetric positive semi-denite damping/absorption matrix, and Dsf describes the uid structure coupling.

It should be noted that for the uid model, the nite element discretization applied to the weak form of the partial differential equation is used. In contrast to this in the model describing the vibration of the structure, in industrial and engineering practice, typically a direct discrete nite element model is employed. This leads to one of the challenges that we will discuss below.

For the displacement vector ud in the different nite elements of the discretized structure, assuming linear material laws, one directly obtains a linear second order system of differential-algebraic equations given by

Msd + Ds ud + Ksud = fe + fp, (4)

where fe is a (discrete) external load and fp is the pressure load. Here Ms is a symmetric positive denite mass matrix, and Ds is a symmetric positive semi-denite damping matrix; both are real. The matrix Ks has the form Ks = K1() + K2 with

real symmetric K1, K2, where K1 is the positive semi-denite stiffness matrix. It is often frequency dependent to model nonlinear material behavior. The matrix K2 models hysteretic damping, that is, damping that is proportional to the displacement (instead of the velocity), but is in phase with velocity [6]. Typically, the matrix K2 is of very small rank.

Although Ms is positive denite in theory, the matrix encountered in practice is highly singular due to the fact that rotational masses are omitted. On the positive side, it is block diagonal with small blocks.

It remains to further discuss the coupling of the uid part and the structure part of the system. The term fp originates from the pressure load Fp =

S p dS. Finite

Page 6 of 18 Mehrmann, Schrder

element modeling/discretization yields fp = DTsf pd and hence,

Msd + Ds ud + Ksud DTsf pd = fe. (5)

If we combine all the equations then we obtain a second order system of differential-algebraic equations

[bracketleftBigg]

Ms 0

DTsf Mf

[bracketrightBigg] [bracketleftbigg]

ud

pd

[bracketrightbigg] + [bracketleftbigg]

Ds 0

0 Df

[bracketrightbigg][bracketleftbigg]

ud

pd

[bracketrightbigg] + [bracketleftbigg]

Ks() DTsf

0 Kf

[bracketrightbigg][bracketleftbigg]

ud pd

[bracketrightbigg] = [bracketleftbigg]

fs

0

[bracketrightbigg]

. (6)

Typically in this model the format of Ms is a factor 1,000-10,000 larger than the format of Mf , the mass and stiffness matrices are highly dependent on the geometry and topology of the car, and the type of nite elements that are used. Since the structure is essentially modeled with a ne (pretty much uniform) mesh, the matrices have dimensions of several millions. The stiffness matrix depends also on the material parameters, while the damping and coupling matrices depend on geometry, topology and material parameters.

2 Results and discussion

2.1 Direct frequency response

Usually one of the rst tasks in the analysis of acoustic elds is to solve the frequency response problem, that is, to simulate the behavior of the acoustic eld under excitations (typically of the structure). The classical approach for the frequency response analysis of linear vibrational systems is to perform a Fourier ansatz

[bracketleftbigg]

ud pd

[bracketrightbigg] = [bracketleftbigg]

u

p

[bracketrightbigg]

et, fs = f ()et,

which then leads to a frequency dependent linear system of equations

2 [bracketleftbigg]

Ms 0 DTsf Mf

[bracketrightbigg] +

Ds 00 Df[bracketrightbigg]

(7)

This linear system, which (for a reasonable structure) has several millions of equations, typically has to be solved for a large frequency range = 0, . . . , 1,000 Hz in

small frequency steps and, depending on the type of excitation, for many right hand sides (load vectors).

Based on the frequency response analysis it is then possible to detect places where the excitation leads to large noise emissions (hot spots), and this approach can be used to improve or even optimize the frequency response behavior within an optimization loop.

+

Ks() DTsf 0 Kf

[bracketrightbigg][parenrightbigg] [bracketleftbigg]

u()

p()

[bracketrightbigg] = [bracketleftbigg]

f ()

0

[bracketrightbigg]

.

Journal of Mathematics in Industry (2011) 1:7 Page 7 of 18

Before we discuss solution methods let us introduce a useful modication of (7). For nonzero frequencies we may multiply the second block row by 1 and introduce the new variable v() = 1 p(). Then we obtain a new linear system

2 [bracketleftbigg]

Ms 00 Mf[bracketrightbigg] +

[bracketleftBigg]

Ds DTsf DTsf Df

[bracketrightBigg]

(8)

which now has complex symmetric coefcients and, hence allows the use of memory-and op-saving structured methods. To simplify notation, we write this as

P ()x() :=

2M + D + K()

x() = b(). (9)

Since this system of equations has to be solved for many right hand sides and over a large frequency range, and all this within an optimization loop, a very efcient solver is necessary. This solver has to be able to recycle information from nearby problems (when stepping through the frequency range and modifying parameters in an optimization) and has to work efciently in a modern multi-processor, multi-core hardware environment. Altogether, this is really a lot to ask for from a linear system solver and makes the use of black box solvers extremely difcult.

2.1.1 State-of-the-art in linear system solvers

Modern techniques for the solution of linear systems in industrial applications are typically a combination of methods, that use the best of each of the various available classes of methods.

The rst class of methods are highly efcient direct solution methods that use Gaussian elimination with partial pivoting combined with other techniques from graph theory to achieve optimal performance, make efcient use of the sparsity structure to avoid ll-in, and save storage. Many of them are even implemented for the use on multiprocessor machines. Well known packages include UMFPACK [7], PARDISO [8], MUMPS [9], WSMP [10], or the HSL collection (for example, HSL_MA78 [11]) to name a few. In view of possibly many right hand sides they are a clear option, despite the fact that for the size and type of problems considered in industrial practice, the storage demands are so high that they typically can only work in an out-of-core setting, that is, the matrix factors are stored on hard disk instead of main memory. Another difculty is that a new factorization has to be computed for every frequency or every modication in the optimization loop. Finally also the bad scaling and the fact that the problems get increasingly ill-conditioned for large frequencies presents a real challenge, because the desired solution accuracy may not be realizable. Some of the packages provide specialized routines for complex symmetric systems, which have the potential to half the storage and computational requirements, but may also increase the complexity of pivoting strategies.

The second class of methods are iterative methods of the Krylov subspace type like the generalized minimal residual method (GMRES) [12], the Bi-conjugate gradient method in its stabilized form (BICGSTAB) [13], or the quasi-minimal residual

+

Ks() 0 0 Kf

[bracketrightbigg][parenrightbigg] [bracketleftbigg]

u()

v()

[bracketrightbigg] = [bracketleftbigg]

f ()

0

[bracketrightbigg]

,

Page 8 of 18 Mehrmann, Schrder

method (QMR) [14]. A variant of the conjugate gradient method for complex symmetric systems like (9) is CSYM [15]. In principle these methods would be very good candidates, since the linear systems are sparse and matrix vector multiplications are relatively cheap. Also, if f and K are independent of , then one can exploit the

shift invariance property of Krylov subspaces [16, 17]. However, often f or K do

depend on and without a good preconditioner the convergence of iterative methods (in particular for large frequencies) is dramatically slow or not realizable.

The third class of potential methods (exploiting the fact that there is a partial differential equation in the background) are multigrid or multilevel methods [1820]. Unfortunately, despite the fact that they are efciently used in the solution of wave propagation problems, such as Helmholtz or Maxwell equations [21], they cannot be used in a simple way in the described acoustic eld problems, since the discrete FEM modeling of the structure does not provide a nice hierarchy of basis functions. If such a hierarchy was available then these methods would lead to very good preconditioners for the Krylov subspace methods. However, with the current state of discrete FEM modeling (being engineering practice for decades) it is extremely difcult to incorporate them into the current code solvers. This is even true for the algebraic multigrid methods like the Ruge-Stben method [22] or hypre [23], due to the difculty of having to solve for a large frequency range, many right hand sides and all this inside an optimization loop. It is needless to say that also the incorporation of these solvers in a multiprocessor, multi-core industrial software tool is another challenge.

The fourth class of methods which have made tremendous progress in the last decades are the adaptive nite element methods, see, for example, [24]. These rene the computational grid according to a priori and a posteriori error estimates and if implemented properly avoid globally ne meshes and therefore the extremely large linear systems and eigenvalue problems in the rst place. However, the construction, analysis, and implementation of such methods for their use inside an industrial package, including the treatment of a full frequency range and many right hand sides, is still in its infancy and requires a major research effort which fortunately is currently addressed in several research projects world-wide.

When our cooperation with SFE started, we essentially made the above analysis of the available methods and realized the following major obstacles.

The problems are badly scaled and get increasingly ill-conditioned when grows; for some parameter constellations the system becomes exactly singular with pos

sibly inconsistent right hand sides;

typically there are many right hand sides; classical off-the-shelf iterative methods do not work well, or fail completely; direct solution methods have to work by storing the factors in out-of-core memory

and cannot easily recycle information from previous frequency or optimization steps;

no multilevel or adaptive grid renement is applicable, the methods must be matrix

based;

the matrices M, K1 are often slightly indenite because of rounding errors during

their creation.

Journal of Mathematics in Industry (2011) 1:7 Page 9 of 18

2.1.2 What did we do?

With all the obstacles of the problem and all the deciencies of the current methods, the development and implementation of an appropriate linear solver for acoustic eld frequency response problems is a major challenge. Furthermore, as in all industrial projects, there was a deadline. In such a circumstance, the only possibility is to compromise between the optimality of a given method for a given problem, the provability of success, and the practical needs in the industrial environment.

Our compromise, since we had to deal with given matrices (as opposed to partial differential equations), was to develop and implement (in the SFE software environment) a preconditioned subspace recycling Krylov subspace method which has the following features:

The initial search space for any frequency is chosen to contain the solutions for the last few frequencies, an approach also used in [25]. For the preconditioner, whose construction and application represent the dominant cost of the algorithm, we had to fulll the requirement that it had to run in a distributive setting, store its results out-of-core, support complex arithmetic, in particular support complex symmetric systems. This is a tall order to ask for but fortunately we were able to base the preconditioner on existing software, using the complex LDLT factorization code of the MUMPS package [9]. The preconditioner was constructed from the symmetric real part of the linear system (9), that is, from

P () := 2 [bracketleftbigg]

For small values of up to ca 60 2 corresponding to a frequency of 60 Hz, typ

ically only 3-6 iteration steps per frequency were necessary and the preconditioner could be kept unchanged for the whole frequency range. But the number of iteration steps per frequency grows substantially the larger gets, so that more and more new preconditioners have to be computed. For large frequencies close to 1,000 Hz, the preconditioned Krylov method tends to be extremely slow or not convergent at all. As a consequence we constructed a hybrid method, where a certain frequency onwards MUMPS is used as a direct solver for the system with the full matrix P (). Furthermore, it was decided to make use of the multi-core environment and to treat several linear systems simultaneously in a distributed fashion. One group of processors solves the systems corresponding to the lower frequencies using the described iterative recycling method. Meanwhile, the other processor groups treat the high frequency systems with a full direct solve for each value of .

In Figure 2 we present computation times of our solver for a model with 551,388 degrees of freedom, solved on a 4 socket machine equipped with Intel Xeon X5670 processors (resulting in 48 virtual CPUs) using 7.9 GB of RAM under Linux 2.6.34. All factorizations were performed out-of-core. It can be observed that computation times scale almost perfectly with the number of systems treated concurrently and that an increase in the number of right hand sides causes only a mild increase in computation time. Moreover, times can be greatly reduced by distributing the solution (LDLT decomposition) of every individual system on multiple processors, although in this case scaling is not perfect.

Ms 0

0 Mf

K1 0

0 Kf

[bracketrightbigg] + [bracketleftbigg]

[bracketrightbigg]

.

Page 10 of 18 Mehrmann, Schrder

Fig. 2 Computation times of the described solver for 1, 2, 4, and 8 concurrent freqencies with different number of processors per frequency (1 or 6) and different number of right hand sides (RHS).

2.1.3 Evaluation of our approach

Since there are commercial packages available, a question that may arise is why we need new numerical methods and solvers in the rst place. From the point of view of the company SFE this is clear, they wanted their own solver in their product, and the solution method including the implementation was satisfactory for their needs. But in many respects the current methods are unsatisfactory from a scientic point of view. The discussed industrial problems present some of the current grant challenges in the eld and the commercial packages are by no means able to solve these problems in a completely satisfactory way. To always get a more-or-less reasonable solution, and to be efcient, they often compromise for accuracy of the results which is certainly inadequate from a scientic point of view. Thus the cooperation with the industrial partner triggered new research questions for the academic world.

It became clear during the project that the concept of recycling in iterative methods, that is, the reuse of information that was already computed for other frequencies or in the course of an optimization procedure is not well-enough understood. As a consequence, motivated by the project with SFE, a new research project [26] in the DFG Research Center MATHEON was started to further investigate the basic mathe

matical principles and to understand how they can be incorporated into new efcient methods [27].

As a second major obstacle, we identied the use of discrete FEM modeling in structural engineering. It would be much easier if adaptive FEM would be usable for the discussed class of problems, and it would also be good to have a grid hierarchy that allows the use of efcient multilevel preconditioners [1820]. Despite the high research activity in this eld this is a major challenge for the described problem classes and almost no analysis or methods are available. We will discuss this in more detail in the section devoted to eigenvalue computation, but again here we see a need for a stronger cooperation between academia and industry in this area of structural engineering, to transfer new ideas that are developed now into the industrial practice.

Journal of Mathematics in Industry (2011) 1:7 Page 11 of 18

2.2 Modal reduction

We have seen in the previous section that due to the ne mesh used for the discrete FEM modeling of the structure, the linear systems have very large dimensions. Typically, however, one is interested in damping the low frequencies associated with the eigenvalues in the neighborhood of 0 (and the imaginary axis) of the complex symmetric matrix function

Q() := 2M + D + K = Q()T , (10) with M, D, K as in (9). The methods discussed below assume that P is quadratic in for the eigenvalue computation, that is, the nonlinear dependency of K on the frequency is ignored. One easily veries that if 0 is an eigenvalue of Q(), that is,

P (0)x = 0, then xT is the corresponding left eigenvector but essentially no other

general properties of the eigenvalue problem are available, in particular, there is unfortunately no immediate variational property for the eigenvalues/eigenfunctions available, as there is for the undamped case [28].

Since one is interested only in part of the spectrum, a natural idea is to identify a space associated with eigenvectors and generalized eigenvectors associated with the important eigenvalues, and to project the problem into this subspace. This is a model reduction approach called modal reduction which, if efciently implemented, can save a huge amount of computing time and storage, despite the fact that it is partially heuristic.

2.2.1 Complex symmetric quadratic eigenvalue problems

In view of the desire to do modal reduction, another task in the cooperation with SFE GmbH was the computation of a small number, say , eigenvalues in a trapezoidal region near zero, see Figure 3(left), and the corresponding subspace S spanned by

the eigenvectors and generalized eigenvectors associated with these eigenvalues, that is, S = span{s1, . . . , s }, where the si form an orthonormal basis of this subspace,

that is, S = [s1, . . . , s ] is an isometry.

Fig. 3 A typical trapezoidal region - within which all eigenvalues are sought - at beginning and after three shifts have been processed.

Page 12 of 18 Mehrmann, Schrder

The projected system would have the form

Q () := 2M + D + K := 2ST MS + ST DS + ST KS (11)

and would still be complex symmetric. In the context of model reduction, the requirements for the reduced model are that the projected system is a good approximation to the large scale system for a large frequency range, and also for a large set of parameter variations. Furthermore, it has to be of small enough dimension, so that classical methods for nonlinear non-smooth optimization can be applied. Again this is a lot to ask, since currently for large scale problems there are really no methods available that guarantee to obtain all the eigenvalues of (10) in a specied region of the complex plane. For small dense problems one could employ the sign function method, [29, 30] but this would require storing full dense inverses of matrices of the given class, which is certainly not possible in the described acoustic eld problems.

The currently used industrial techniques typically solve a simplied problem, for example, the eigenvalues and associated eigenvectors in the desired region for the undamped/uncoupled problem are used for the projection or an algebraic multilevel substructuring method (AMLS) is used that exploits the structure of the matrices as they arise form the discrete FEM [31]. All these techniques are partially heuristic, since there is no guarantee that all the desired eigenvalues are captured or that the eigenvalues from the projected system are close to those of the original physical system or of the FEM discretized system, that is, in general, currently no error bounds are available. To generate such error bounds is a major challenge for the academic community and research in this direction has started in the project [32] with rst results in [33, 34].

Furthermore, the numerical solution of the large scale nonlinear eigenvalue problem (10) itself also presents a mathematical challenge, for several reasons.

First of all the mass matrix is singular, so the problem has eigenvalues at which

often cause convergence problems for the iterative methods, second of all, there are also several (typically six) eigenvalues at 0, corresponding to the six free degrees of motion, three each for translation and rotation, of the whole structure. In industrial practice the eigensolver is furthermore also often used for model verication, by checking whether there are exactly 6 eigenvalues at 0. Then, if the model is awed, the singularity may be even higher.

The major challenge, however, is that we want the eigenvalues near 0 and it is well known that classical iterative methods for large sparse eigenvalue problems, like the implicitly restarted Arnoldi method [35], typically converge fast only to the eigenvalues at the periphery of the spectrum [36]. Thus, either a shift-and-invert technique that transforms the problem and maps the desired eigenvalue to the periphery is necessary or other techniques like Newtons method [37, 38] or the Jacobi-Davidson method [39, 40] have to be used. All these methods require the solution of linear systems of the form (

2M +

D + K)x = b with a given shift-point

. But this is the problem that we wanted to avoid in the rst place and hence a vicious circle is closed.

Journal of Mathematics in Industry (2011) 1:7 Page 13 of 18

2.2.2 What did we do?

Having assessed the various options and their potential advantages and drawbacks, we decided to develop a new method based on the following concepts and to implement it into the SFE environment.

First of all, since the methods that work directly on the quadratic eigenvalue problems, for example, [28, 41, 42], are not yet as mature as those for linear eigenvalue problems, the problem is linearized by introducing a new variable x and turning (10) from a quadratic into a linear eigenvalue problem.

Since there is no method available that is guaranteed to nd all eigenvalues in a region at feasible costs, we had to revert to a partially heuristic approach. We are looking for eigenvalues in the vicinity of target points which are scattered inside the region of interest. To nd eigenvalues close to a particular target , also called a shift, we use the shift-linearize-invert-ansatz [42] that requires computation of the largest eigenvalues of the matrix

A :=

This is done by the block Arnoldi method [43] which requires the application of the matrix A to given blocks of vectors during the generation of the Krylov subspace.

This involves a solution of a linear system with the complex symmetric matrix Q(). We again use the direct solver MUMPS to compute a complex LDLT factorization.

The block Arnoldi method searches for eigenvectors of A within the block Krylov

subspace

Km(A, B) := span

B, AB, A2B, . . . , Am1B

, (12) which is known to contain, for increasing m, increasingly good approximations to eigenvectors corresponding to eigenvalues of A of maximum absolute value - which

correspond to the eigenvalues of Q() closest to . In the classical Arnoldi method, the starting block B is a vector. Otherwise, when B consists of several, say nb, columns, the resulting method is the block Arnoldi method. The matrix B can be chosen in an (almost) arbitrary way, but if eigenvector or subspace approximations are known, the use of these will speed up convergence. Among the advantages of using the block method are that clusters of eigenvalues are handled better [44]. Furthermore, all vector-vector operations become matrix-matrix operations which can be implemented much more efciently by use of BLAS level 3 routines [45] and the matrix-vector multiplication with the large matrix A can be carried out with a block

of vectors. Furthermore, the heuristic part of the algorithm is more reliable in the block case, see below.

In the actual implementation of the block Arnoldi algorithm the terms AiB are

never explicitly formed. Instead an orthonormal basis Qm = [V1, V2, . . . , Vm] of Km(A, B) is generated by setting V1 to be an orthonormal basis of span(B) fol

lowed by iteratively forming Vi+1 by orthonormalization of AVi against V1, . . . , Vi.

Classically, a variant of the Gram-Schmidt method is used for this task [35], but any other orthonormalization procedure may be employed. We use block Householder

Q()1 0 0 I

[bracketrightbigg][bracketleftbigg]2M

+ D M I 0

[bracketrightbigg]

.

Page 14 of 18 Mehrmann, Schrder

reections [46, 47] that are very stable and attain high efciency by using BLAS-3 operations. Collecting all the Vi and the orthonormalization coefcients Hi,j results in the well known Arnoldi relation AQm = QmHm + Vi+1Hm+1,mETm with an

(m m)-block Hessenberg matrix Hm = (Hi,j )mi,j=1. The eigenvalues of Hm, called

Ritz values, are used as approximations of eigenvalues of A. Likewise, the eigenvec

tors of Hm, multiplied with Qm, are used as approximate eigenvectors of A.

A drawback of the block Arnoldi method is the necessity to store the basis Qm which grows for increasing m. A typical way out of this problem is to restart the algorithm at some point. Instead of restarting with a single new starting block it is possible to restart with a whole Arnoldi relation. For the vector Arnoldi method, two such implicit restarting schemes are common, one using a lter polynomial [35] and one using a reordering of the Schur form of Hm [48]. We are using the latter approach as the rst one does not elegantly generalize to the block case [46].

The heuristic part of our method is that we assume that the block Arnoldi method really nds all eigenvalues in a circle around the shift . While quite often the eigenvalues do indeed appear and converge in the order of the distance to the shift, it is not rare that one or a group of eigenvalues converge slower than other farer away eigenvalues. However, in these cases usually the missed eigenvalues are present as (yet) unconverged Ritz values. Therefore, we use the unconverged Ritz value that is closest to the shift as the radius of a circle that is trusted to contain no missed eigenvalues.

In the implementation, we repeatedly run the block Arnoldi method for different shifts, possibly several at once in a distributed setting. Figure 3(right) shows the situation after a few iterations. Large parts of the trapezoidal region are covered, leaving only some small remaining regions to be searched. New shifts are placed inside the largest such white regions, until the whole trapezoidal region of interest is covered by trusted circles.

Of course, with such a covering approach an eigenpair could be computed by more than one Arnoldi run for different shifts. For that reason the freshly discovered eigenpairs have to be checked for being copies of already previously found pairs. To achieve this we consider a new eigenpair (, x) being a copy, if [xT , xT ]T is

almost linearly dependent to the span of the vectors [xTold, oldxTold]T corresponding

to every known eigenvalue old sufciently close to .

We applied the solver to a car model, discretized by a regular mesh of 35 mm leading to 219,432 degrees of freedom. We were looking for eigenvalues within the triangle bordered by the lines Im() > 20 Re(), Im() > 20 Re(), and Im() < f 2,

for f {50, 100, 150, 200, 250}. Table 1 lists the number of eigenvalues within these

triangles, the number of used shifts, that is, used matrix factorizations, and the number of overall Arnoldi iterations, that is, matrix-block products. The test was performed on a PC with an Intel Core2 Duo E6850 CPU clocked at 3.0 GHz, with 4 GB RAM. One shift was addressed at a time using one processor. The block size was 5.

2.2.3 Evaluation of our approach

Our approach certainly does not use many novel ideas, but instead builds on mature and proven concepts. The method can run in a distributed setting processing several shifts at once, each one running in several processes. The matrix factorizations can

Journal of Mathematics in Industry (2011) 1:7 Page 15 of 18

Table 1 Number of required shifts and iterations to nd all eigenvalues between 0 and {50, . . . , 250} Hz.

0 Hz to: 50 Hz 100 Hz 150 Hz 200 Hz 250 Hz

Number of found eigenvalues 20 52 129 217 346

Number of shifts 1 3 3 5 7

Number of iterations 39 139 136 294 437

be kept out-of-core. Moreover, since the basis vectors are explicitly orthogonalized by a stable Householder scheme, the generated block Arnoldi basis is orthonormal to working precision. This is in sharp contrast to the basis generated by the non-symmetric Lanczos method which typically looses linear independence after enough iterations and leads to so-called spurious eigenvalues [49, 50].

In our experiments the heuristic choice of the radii of trusted circles worked very well. In thousands of test cases with problems from SFE as well as randomly generated problems, it happened only once that this approach missed an eigenvalue. This happened with a block size of one, that is, with the vector Arnoldi method. The method worked ne for this case if a block size of at lest two was used. In our implementation the default block size is 8.

Unfortunately, the improved robustness of the block-Arnoldi method compared with the standard Arnoldi method comes at the price of increased memory requirements to store the basis. This downside is somewhat mitigated, however, by the use of restarts.

As another downside of our approach, so far only quadratic eigenvalue problems can be solved, for truly nonlinear problems with dependence in K the methods is not directly applicable. On the other hand, the method can be applied to uid or structure subsystems separately, or to the complete coupled system, and it tolerates a singular mass matrix.

3 Conclusions. The two-way street of industrial cooperation

We discussed the modeling, simulation and (at least the vision) of the optimization of acoustic elds. As an example of our industrial cooperation we studied frequency domain computations and modal reduction methods for the acoustic eld inside a car, as well as the challenges of the resulting linear systems and eigenvalue problems. One lesson from the project was that the methods to be used in industrial practice cannot be constructed and implemented using textbook approaches. For instance, very often tricks have to be employed that increase the efciency of the computation, but that lack a full mathematical understanding. This can lead to misunderstandings between engineers, programmers and mathematicians. A stronger communication and cooperation between these groups is necessary to address the described challenges. If this is achieved then all sides benet from a cooperation. The work with SFE GmbH on the frequency response problem for interior acoustic eld computation started out as a clear transfer project (a one way street), with the idea to transport the knowledge and know-how about current linear system and eigenvalue solver technologies available in the academic environment into an industrial software environment.

Page 16 of 18 Mehrmann, Schrder

But as we have discussed, already early on in the cooperation a lot of new research topics appeared that could not be treated within the current project (two-way-street). Examples include recycling methods for a sequence of slowly changing linear systems, updating of preconditioners, guaranteed location of all eigenvalues inside a region of the complex plane, to name a few.

Competing interests

Our research was supported by Deutsche Forschungsgemeinschaft, via the DFG Research Center MATH-

EON, Mathematics for Key Technologies, in Berlin and by SFE GmbH, Voltastr. 5, 13355 Berlin.

Authors contributions

All authors participated in designing and implementing the solvers as well as in writing the manuscript. All authors read and approved the nal manuscript.

Acknowledgements We thank H. Zimmer and A. Hilliges from our industrial partner SFE for the permission to use the gures obtained in the cooperation, and we thank Tobias Brll, Elena Teidelt, Matthias Miltenberger for their help in the implementation. Moreover, we thank two anonymous referees for their valuable comments on a preliminary version of this paper.

References

1. ANSYS, Inc.: Programmers manual for ANSYS. http://www1.ansys.com/customer/content/documentation/110/ansys/aprog110.pdf

Web End =http://www1.ansys.com/customer/content/ http://www1.ansys.com/customer/content/documentation/110/ansys/aprog110.pdf

Web End =documentation/110/ansys/aprog110.pdf (2007)

2. MSC. Software Corporation: MSC.Nastran 2005 Installation and Operations Guide. http://www.mscsoftware.com/Products/CAE-Tools/MSC-Nastran.aspx

Web End =http://www. http://www.mscsoftware.com/Products/CAE-Tools/MSC-Nastran.aspx

Web End =mscsoftware.com/Products/CAE-Tools/MSC-Nastran.aspx (2004)

3. Antoulas, A.C.: Approximation of large-scale dynamical systems. In: Advances in Design and Control, vol. 6. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA (2005)

4. Benner, P., Mehrmann, V., Sorensen, D. (eds.): Dimension Reduction of Large-Scale Systems. LNSCE, vol. 45. Springer, Heidelberg (2005)

5. Schmidt, F., Deuhard, P.: Discrete transparent boundary conditions for the numerical solution of Fresnels equation. Comput. Math. Appl. 29, 5376 (1995)

6. Bishop, R.: The treatment of damping forces in vibration theory. Jl R. Aeronaut. Soc. 59, 738 (1955)7. Davis, T.A.: Algorithm 832: UMFPACK V4.3 - an unsymmetric-pattern multifrontal method. ACM Trans. Math. Softw. 30(2), 196199 (2004)

8. Schenk, O., Grtner, K.: Solving unsymmetric sparse systems of linear equations with PARDISO. Future Gener. Comput. Syst. 20(3), 475487 (2004). http://www.sciencedirect.com/science/article/B6V06-49NXY7J-F/2/e8260e5d8f19639019cddea4776c024c

Web End =http://www.sciencedirect.com/ science/article/B6V06-49NXY7J-F/2/e8260e5d8f19639019cddea4776c024c

9. MUMPS: MUltifrontal Massively Parallel Solver: Users Guide. http://graal.ens-lyon.fr/MUMPS

Web End =http://graal.ens-lyon.fr/MUMPS (2009)

10. Gupta, A.: WSMP: Watson Sparse Matrix Package, Part II. Direct solution of general sparse systems. IBM research report, IBM T. J. Watson Research Center. http://www.alphaworks.ibm.com/tech/wsmp

Web End =http://www.alphaworks.ibm.com/tech/wsmp (2000)

11. Reid, J.K., Scott, J.A.: An efcient out-of-core multifrontal solver for large-scale unsymmetric element problems. Technical Report RAL-TR-2007-014, SFTC Rutherford Appleton Laboratory (2007)

12. Saad, Y., Schultz, M.H.: GMRES: a generalized minimal residual algorithm for solving nonsymmetric linear systems. SIAM J. Sci. Stat. Comput. 7(3), 856869 (1986)

13. van der Vorst, HA: Bi-CGSTAB: a fast and smoothly converging variant of Bi-CG for the solution of nonsymmetric linear systems. SIAM J. Sci. Stat. Comput. 13(2), 631644 (1992)

Journal of Mathematics in Industry (2011) 1:7 Page 17 of 18

14. Freund, R.W., Nachtigal, N.M.: QMR: a quasi-minimal residual method for non-Hermitian linear systems. Numer. Math. 60(3), 315339 (1991)

15. Bunse-Gerstner, A., Stver, R.: On a conjugate gradient-type method for solving complex symmetric linear systems. Linear Algebra Appl. 287(1-3), 105123 (1999)

16. Gu, G.D., Simoncini, V.: Numerical solution of parameter-dependent linear systems. Numer. Linear Algebra Appl. 12(9), 923940 (2005)

17. Simoncini, V., Perotti, F.: On the numerical solution of (2A + B + C)x = b and application to

structural dynamics. SIAM J. Sci. Comput. 23(6), 18751897 (2002) (electronic)18. Hackbusch, W.: Multigrid Methods and Applications. Springer Series in Computational Mathematics, vol. 4. Springer-Verlag, Berlin (1985)

19. McCormick, S.F.: Multilevel Adaptive Methods for Partial Differential Equations. Frontiers in Applied Mathematics, vol. 6. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA (1989)

20. McCormick, S.F.: Multilevel Projection Methods for Partial Differential Equations. CBMS-NSF Regional Conference Series in Applied Mathematics, vol. 62. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA (1992)

21. Ernst, O.G.: Fast numerical solution of exterior Helmholtz problems with radiation boundary condition by imbedding. Ph.D. thesis, Stanford University (1994)

22. Ruge, J.W., Stben, K.: Algebraic multigrid. In: Multigrid Methods. Frontiers Appl. Math., vol. 3, pp. 73130. SIAM, Philadelphia, PA (1987)

23. Falgout, R.D., Jones, J.E., Yang, U.M.: The design and implementation of hypre, a library of parallel high performance preconditioners. In: Numerical Solution of Partial Differential Equations on Parallel Computers. Lect. Notes Comput. Sci. Eng., vol. 51, pp. 267294. Springer, Berlin (2006). doi:10.1007/3-540-31619-1_8

24. Bangerth, W., Rannacher, R.: Adaptive Finite Element Methods for Differential Equations. ETH Lecture Series. Birkhuser, Basel, Ch. (2003)

25. Fischer, P.F.: Projection techniques for iterative solution of Ax = b with successive right-hand sides.

Comput. Methods Appl. Mech. Eng. 163(1-4), 193204 (1998)26. Numerical methods for large-scale parameter-dependent systems. Project C29, DFG Research Center Matheon, Berlin. http://www.matheon.de/research/show_project.asp?id=172

Web End =http://www.matheon.de/research/show_project.asp?id=172

27. Gaul, A., Gutknecht, M.H., Liesen, J., Nabben, R.: Deated and augmented Krylov subspace methods: Basic facts and a breakdown-free deated MINRES. Preprint 749, DFG Research Center Matheon, Mathematics for key technologies, Berlin (2011)

28. Mehrmann, V., Voss, H.: Nonlinear eigenvalue problems: a challenge for modern eigenvalue methods. Mitt. der Ges. f. Angewandte Mathematik und Mechanik 27, 121151 (2005)

29. Bai, Z., Demmel, J.: Using the matrix sign function to compute invariant subspaces. SIAM J. Matrix Anal. Appl. 19, 205225 (1998). doi:10.1137/S0895479896297719 (electronic)

30. Higham, N.J.: Functions of matrices. Theory and Computation. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA (2008)

31. Bennighof, J.K.: Adaptive multi-level substructuring method for acoustic radiation and scattering from complex structures. In: Kalinowski, A.J. (ed.) Computational Methods for Fluid/Structure Interaction, vol. 178, pp. 2538. AMSE, New York (1993)

32. Adaptive solution of parametric eigenvalue problems for partial differential equations. Project C22, DFG research center Matheon, Berlin. http://www.matheon.de/research/show_project.asp?id=126

Web End =http://www.matheon.de/research/show_project.asp?id=126 ]

33. Carstensen, C., Gedicke, J., Mehrmann, V., Miedlar, A.: An adaptive homotopy approach for nonselfadjoint eigenvalue problems. Preprint 718, DFG Research Center Matheon, Mathematics for key technologies, Berlin (2010)

34. Miedlar, A.: Inexact adaptive nite element methods for elliptic PDE eigenvalue problems. Ph.D. thesis, Technische Universitt Berlin, Inst. f. Mathematik (2011)

35. Lehoucq, R.B., Sorensen, D.C., Yang, C.: ARPACK Users Guide. Solution of Large-Scale Eigenvalue Problems with Implicitly Restarted Arnoldi Methods. Software, Environments, and Tools, vol.6. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA (1998)36. Saad, Y.: Numerical Methods for Large Eigenvalue Problems. Algorithms and Architectures for Advanced Scientic Computing. Manchester University Press, Manchester (1992)

37. Schwetlick, H., Schreiber, K.: A primal-dual Jacobi-Davidson-like method for nonlinear eigenvalue problems. Internal Report ZIH-IR-0613, Techn. Univ. Dresden, Zentrum fr Informationsdienste und Hochleistungsrechnen (2006)

38. Schwetlick, H., Schreiber, K.: Nonlinear Rayleigh functionals. Linear Algebra and Its Applications. In Press, Corrected Proof (2010). doi:http://dx.doi.org/10.1016/j.laa.2010.06.048

Web End =10.1016/j.laa.2010.06.048

Page 18 of 18 Mehrmann, Schrder

39. Sleijpen, G.L.G., Booten, A.G.L., Fokkema, D.R., Van der Vorst, H.A.: Jacobi-Davidson type methods for generalized eigenproblems and polynomial eigenproblems. BIT Numer. Math. 36(3), 595633 (1996)

40. Sleijpen, G.L.G., Van der Vorst, HA: A Jacobi-Davidson iteration method for linear eigenvalue problems. SIAM J. Matrix Anal. Appl. 17(2), 401425 (1996)

41. Beyn, W.: An integral method for solving nonlinear eigenvalue problems. ArXiv e-prints (2010)42. Tisseur, F., Meerbergen, K.: The quadratic eigenvalue problem. SIAM Rev. 43(2), 235286 (2001)43. Arnoldi, W.E.: The principle of minimized iteration in the solution of the matrix eigenvalue problem.Q. Appl. Math. 9, 1729 (1951)44. Lehoucq, R., Maschhoff, K.: Implementation of an implicitly restarted block Arnoldi method. Preprint MCS-P6490297, Argonne National Laboratory, Argonne, IL (1997)

45. Dongarra, J.J., Croz, J.D., Hammarling, S., Duff, I.: A set of level 3 basic linear algebra subprograms. ACM Trans. Math. Softw. 16, 117 (1990). doi:10.1145/77626.79170

46. Baglama, J.: Augmented block Householder Arnoldi method. Linear Algebra Appl. 429(10), 2315 2334 (2008)

47. Schreiber, R., Van Loan, C.: A storage-efcient WY representation for products of Householder transformations. SIAM J. Sci. Stat. Comput. 10, 5357 (1989)

48. Stewart, G.W.: A Krylov-Schur algorithm for large eigenproblems. SIAM J. Matrix Anal. Appl. 23(3), 601614 (2001/02)

49. Cullum, J., Willoughby, R.A.: Lanczos and the computation in specied intervals of the spectrum of large, sparse real symmetric matrices. In: Sparse Matrix Proceedings 1978 (Sympos. Sparse Matrix Comput., Knoxville, Tenn., 1978), pp. 220255. SIAM, Philadelphia, Pa (1979)

50. Freund, R.W., Gutknecht, M.H., Nachtigal, N.M.: An implementation of the look-ahead Lanczos algorithm for non-Hermitian matrices. SIAM J. Sci. Comput. 14, 137158 (1993)