A greedy reduced basis scheme for multifrequency solution of structural acoustic systems

The solution of frequency dependent linear systems arising from the discretization of vibro‐acoustic problems requires a significant computational effort in the case of rapidly varying responses. In this paper, we review the use of a greedy reduced basis scheme for the efficient solution in a frequency range. The reduced basis is spanned by responses of the system at certain frequencies that are chosen iteratively based on the response that is currently worst approximated in each step. The approximations at intermediate frequencies as well as the a posteriori estimations of associated errors are computed using a least squares solver. The proposed scheme is applied to the solution of an interior acoustic problem with boundary element method (BEM) and to the solution of coupled structural acoustic problems with finite element method and BEM. The computational times are compared to those of a conventional frequencywise strategy. The results illustrate the efficiency of the method.


INTRODUCTION
Solving a sequence of linear systems of equations accounts for a significant computational effort that is required for the simulation of time-harmonic vibro-acoustic problems, particularly in the case of rapidly varying responses. Modal superposition is a popular numerical tool to address this issue and often used to avoid solutions for each individual frequency of interest. The respective eigenvalue problem needs to be solved only once and the responses at certain frequencies as well as the responses for different excitations are simply obtained by matrix vector multiplications. Especially in the case of both interior acoustic problems 1 and coupled structural acoustic problems in bounded domains, 2 modal superposition techniques based on the finite element method (FEM) 3 are well established. The respective coefficient matrices are independent of the frequency, and consequently, standard eigenvalue solvers 4 can be applied in a straightforward manner. However, in the case of exterior acoustic problems, the situation is different due to the representation of the unbounded acoustic domain. Three families of methods are available for addressing this issue: FEM with special boundary conditions, 5 the infinite element method (IFEM), 6 and the boundary element method (BEM). 7,8 The use of conjugated Astley-Leis IFEM results in frequency independent coefficient matrices giving rise to the formulation of a generalized eigenvalue problem for determining acoustic radiation modes and normal modes. 9,10 However, making informed choice of the contributing modes for a superposition is still an active field of research. 11 While BEM is particularly advantageous due to the implicit satisfaction of the far field radiation condition, its major drawback stems from the frequency dependence of the coefficient matrices. Special considerations such as the multiple reciprocity method 12 and the contour integral method 13,14 are usually employed to obtain acoustic modes with BEM. In the case of weakly coupled structural acoustic interaction, the superposition of isolated structural modes is sufficient to obtain a reduced model of the underlying vibro-acoustic problem. However, for lightweight structures submerged in heavy fluids, the in-vacuo modes do no longer form the modal basis. Hence, the modal analysis has to be performed including the effect of acoustic loading. When using BEM for the discretization of the fluid, this involves the solution of a nonlinear eigenvalue problem. Peters et al have proposed a frequency approximation of the boundary element impedance matrix 15,16 to overcome this issue, and more recently, the nonlinear structural acoustic eigenvalue problem has been solved using contour integration. 17,18 Despite the abovementioned advances in numerical modal analysis, a frequencywise response analysis is still the most popular procedure for the calculation of vibro-acoustic responses in a frequency range. Moreover, implementations of eigenvalue solvers such as the contour integral method nevertheless require system evaluations at a considerable number of discrete frequency points. Therefore, in the following, we contribute to the acceleration of vibro-acoustic frequency response analysis as an alternative approach to the abovementioned modal superposition techniques.
During the last decades, many researchers have made effort to accelerate the frequency sweep of acoustic as well as structural acoustic problems. Moreover, efficient procedures for the assembly of boundary element matrices in a frequency range have been developed. [19][20][21][22] Approximating the acoustic response itself can also be efficient in cases where only a small partition of the solution is of interest. 23,24 Furthermore, the high numerical complexity associated to fully populated boundary element matrices led to the development of several fast algorithms 25,26 that have also been combined with multifrequency strategies. 27 In addition, a number of publications [28][29][30] are concerned with efficient schemes for the actual solution of linear systems in a frequency range. In particular, projection-based model order reduction (MOR) techniques have been proposed to reduce the computational effort. They have been extensively reviewed for vibro-acoustic analyses using FEM. 28 In the context of coupled FEM-BEM analyses, MOR has been pioneered by the work of Peters et al, 31 in which the structural finite element matrices are reduced by means of Krylov subspace MOR, whereas the frequency-dependent boundary element matrices are left unchanged. More recently, resolvent sampling-based MOR has been proposed for the fully coupled FEM-BEM system. 30 In this method, the projection matrices are spanned by the response at some sampling points inside a frequency interval. Typically, such frequency sampling is done based on the Chebyshev nodes in an interval 15,30,32 without a priori knowledge of the solution and the occurrence of resonances. To the best of our knowledge, there is no strategy available allowing for an optimal choice of frequency samples that goes in hand with an estimation of the error introduced by the projective MOR in the context of vibro-acoustic analyses.
In this paper, we contribute to this subject by proposing a greedy reduced basis scheme for the solution of purely acoustic and coupled structural acoustic problems. In the context of a multifrequency analysis, a reduced basis method 33 expresses the (approximate) solution at each frequency point as a linear combination of a few basis vectors. In our case, these basis vectors are simply the solutions at suitably chosen frequency points within the frequency range of interest. When these basis vectors are chosen iteratively, based on the solution that is worst approximated by the current reduced basis in each step, this procedure is known as a greedy algorithm. 34,35 Section 2 proceeds with a brief review of the coupled FEM-BEM formulation followed by a detailed description of the greedy reduced basis scheme and its implementation in Section 3. The method is then applied to the solution of an interior acoustic problem and to the solution of two coupled structural acoustic problems in Section 4. In all three examples, the greedy algorithm automatically constructs the basis vectors at or near resonances, thus leading to small errors in a few iterations. Section 5 concludes with a discussion on the applicability as well as the current limitations of the method. Generally speaking, the greedy reduced basis scheme proves to be particularly efficient, when the responses at a large number of frequency points are contained in a small subspace. Consequently, the common practice among engineers of oversampling the frequency range has less impact on the overall computational effort than when using a conventional frequencywise strategy. On the other hand, when the responses at different frequencies are not or hardly related to each other (eg, in the case of frequency dependent forcing vectors), conventional frequencywise strategies or modal superposition techniques could be more efficient.

COUPLED FEM-BEM FORMULATION FOR STRUCTURAL ACOUSTIC INTERACTION
We consider fully coupled structural acoustic interaction problems using FEM 36 and BEM 8 for discretizing the equations of linear time-harmonic elastodynamics and acoustics, respectively. The resulting linear systems of equations for the structural and for the acoustic subdomain read The column vectors u ∈ ℂ n s ×1 and p ∈ ℂ n f ×1 contain the unknown displacement and sound pressure degrees of freedom at the nodes. The stiffness and mass matrices of the structure are denoted as K ∈ ℝ n s ×n s and M ∈ ℝ n s ×n s . Structural damping is not considered in this work. The boundary element matrices H( ) ∈ ℂ n f ×n f and G( ) ∈ ℂ n f ×n f are obtained by a collocation discretization of the Kirchhoff-Helmholtz integral equation, relating the structural particle velocity v s ∈ ℂ n s ×1 to the sound pressure. These matrices are implicitly dependent on the angular frequency = 2 f, where f is the frequency in Hz. The structure is excited by nodal forces f s ∈ ℂ n s ×1 . Equations (1) and (2) are coupled on the sound radiating boundary. There, the structure is subject to normal tractions due to the acoustic sound pressure, and the particle velocity in (2) is equal to the time derivative of the normal displacement on the boundary. The coupling conditions can be expressed as where C sf ∈ ℝ n s ×n f and C fs ∈ ℝ n f ×n s are the mesh coupling matrices obtained by a Galerkin projection. 37 The force vector f f can be interpreted as the acoustic loading on the structural nodes and i denotes the imaginary unit. Finally, the global system of equations containing the coupling conditions emerges as When only the structural response is of interest, (4) can be reformulated by forming the Schur complement and thereby omitting the pressure degrees of freedom, 16 ie, Note that, often, the dimension of the acoustic subproblem is at least an order of magnitude smaller than the structural subproblem. Thus, when facing moderately large problems, the computation of H( ) −1 G( ) is relatively inexpensive compared to the solution of (5). However, the reformulation of (4) into (5) deteriorates the sparsity pattern of the structural subsystem and hence precludes the solution of large-scale problems. In these cases, the solution of (4) using a preconditioned iterative scheme is more practical. Alternatively, (5) could be solved using a direct method in conjunction with MOR of the structural subsystem. 31 In the numerical examples in Section 4, we will consider both types of structural acoustic equation (4) and (5). Many vibro-acoustic applications require a solution of (4) or (5) at a range of frequencies 1 , … , m . Usually, each of the m linear systems are solved successively in an independent manner. In the proceeding sections, we propose an alternative approach based on a greedy algorithm.

A GREEDY ALGORITHM FOR MULTIFREQUENCY ANALYSIS
First introduced in the context of combinatorial optimization problems, the general idea behind greedy algorithms is to repeatedly make the choice that is currently optimal. A prominent example for which locally optimal choices (the greedy choices) yield the globally optimal solution is the so-called minimum spanning tree problem. It is the task of connecting a number of vertices such that there exist routes among all vertices while minimizing the total length of the connecting edges. It can simply be proven that this problem can be solved by the following greedy strategy 38 : "In each step, add the shortest edge between two vertices that does not form a closed cycle with the current set of edges." Beside the field of combinatorial optimization, greedy algorithms have nowadays become a popular tool for the construction of reduced bases in the context of parameterized partial differential equations. 39 As such, they represent an alternative to proper orthogonal decomposition (POD) methods, in which possibly large numbers of high-fidelity systems need to be evaluated in order to deduce adequate POD bases. Using greedy algorithms, the number of evaluations of the high-fidelity system is minimized by iteratively adding new basis vectors to the reduced basis. 40,41 However, they rely on an a posteriori error estimation for the whole parameter space in each iteration in order to optimally choose the next parameter sample. Consequently, their success is heavily dependent on the computational cost that is required for the solution of this optimization problem and the associated error estimation. In the following, we propose to extend the application of greedy algorithms to the solution of vibro-acoustic problems in a frequency range. A detailed description of the algorithm is provided including a cheap a posteriori error estimation giving rise to fast frequency range solutions.
For the sake of readability, the frequency dependent linear systems (4) and (5) are replaced by A( )x( ) = b( ) with the frequency points of interest ∈  = { 1 , … , m }. Note that here, we consider the general case with a frequency dependent right-hand side b( ). In each iteration j of the greedy algorithm, a set of frequency samples  ⊆  is given as The corresponding reduced basis X j is simply the concatenation of the solutions at these frequency samples, ie, The current approximation at an arbitrary frequency point ∈  is then expressed by with y( ) ∈ ℂ ×1 , where y( ) is the solution of the least squares problem The main idea behind the greedy algorithm is to choose the next sample (j+1) there, where the approximation (8) yields the largest relative residual. Therefore, holds and the next basis vector is calculated by solving the system The procedure is repeated until a convergence criterion based on the relative residual is fulfilled for all ∈ .
Regarding the computational efficiency of the method, it is crucial that a sufficiently accurate approximation of the solution in the whole frequency range is achievable by a linear combination of only a few q ≪ m solutions at sample frequencies, ie, with a small number of iterations q. From an algebraic point of view, this means that the frequency range solution X = [x( 1 ), … , x( m )] ∈ ℂ n×m admits a low rank approximation. More precisely, the singular values of X should exhibit exponential decay, which has been proven for an analytical frequency dependency of A( ) and b( ) 42 and thus holds for the herein presented case. From a vibro-acoustic point of view, the admissibility of a low rank approximation can be understood in that the spatial distribution of the response exhibits a certain degree of regularity with respect to the frequency. Essentially, the well-known modal reduction technique is based on similar principles of superposition, but instead of superposing modes of the system, here, we use a linear combination of certain responses. In the context of modal superposition techniques, the contribution of the individual modes is expressed by the modal amplification factors. Similarly, using the herein proposed greedy algorithm, the iteratively chosen responses of the system are weighted by the scaling factors y. Although modal quantities are not determined in the present method, we will later see that the frequencies close to resonances are intrinsically chosen as samples over the course of the iterations.
Beside the solution of q linear systems (recall that q corresponds to the total number of iterations), the main computational effort of the greedy algorithm stems from the minimization of (9) using a least squares solver along with an a posteriori error estimation. For both, the built-in MATLAB function lsqlin provides a straightforward implementation. While an error estimation is performed in each iteration for all the yet unconverged solutions and hence necessitate (qm) calls of the least squares solver, the evaluation of a single least squares problem is rather inexpensive since A( )X j has only a small number of columns. The numerical examples will show that the algorithm results in a significantly smaller computational effort than solving m linear systems independently.
The repeated evaluation of the least squares problems, however, requires preassembly and storage of all m possibly fully populated matrices in the main memory. The latter is the major drawback of the proposed method and currently limits its application to small-sized problems. There are two remedies in this regard. For medium-sized problems, the use of data sparse formats such as hierarchical matrices 26 resolves the memory issues. For large-scale problems, the proposed algorithm could be parallelized by solving the least squares problem (9) in a distributed memory environment and solving (11) on a shared memory environment. Although we expect that only a moderate amount of communication is required, systematic benchmarking is still needed to verify the applicability of the proposed scheme to large-scale multifrequency problems. However, this is not within the scope of the present work but planned for future research.
The greedy reduced basis scheme is summarized in Algorithm 1. In order to avoid least squares minimizations at frequency points with already converged solutions, a set of frequencies  sol is defined containing the frequencies at which the solution has been explicitly determined (line 23) as well the frequencies at which the least squares approximation already yields sufficient accuracy (line 15). Furthermore, the matrix products A( i )X j for setting up the least squares minimizations in line 13 are explicitly stored in the main memory to avoid recomputation in each iteration. Only the new column A( i )x( (j) ) is added in each iteration. Since this matrix vector product needs to be evaluated (qm) times, it accounts for the main computational effort beside the explicit solutions of q linear systems and the actual solutions of the (qm) least squares problems.

NUMERICAL EXAMPLES
The proceeding sections verify the efficiency of Algorithm 1 in the context of multifrequency analysis.

Plane sound wave in a closed rigid duct
As a first numerical example, we consider a purely interior acoustic problem. This problem is well suited to demonstrate the convergence behavior of the greedy algorithm in the context of multifrequency analysis. A plane sound wave is excited by a harmonic particle velocity at the left end of a three-dimensional closed and air-filled duct of length l = 3.4 m and square cross section w = 0.2 m. The wave is fully reflected at the acoustically rigid right end. The speed of sound is c = 340 m∕s. The system is free of dissipation, and hence, resonances occur at the integer multiples of 50 Hz. This problem is often used in the computational acoustics community for benchmarking purposes 43 and extensive studies on associated discretization errors with respect to mesh sizes and element types are available in the literature. 44,45 The acoustic field is discretized using a uniform mesh of 1120 quadrilateral boundary elements with bilinear discontinuous sound pressure interpolation yielding 4480 degrees of freedom. We are interested in solving the resultant linear system (2) in the frequency range from 40 to 210 Hz with frequency steps of Δf = 1 Hz, ie, m = 171. Resonances are expected at 50, 100, 150, and 200 Hz. This problem is relatively well conditioned, ie, cond A( i ) ≈ 10 2 .
First, a standard generalized minimal residual (GMRes) algorithm 46 without restarts is used for the purpose of a comparison in terms of the computational time. The tolerance for the relative residual is set to gmres = 10 −5 . All m = 171 linear systems have been successively solved requiring a total of 4463 iterations and a wall clock time of 96.44 s. This corresponds to an average of 0.564 s for the solution at each frequency point.
Using the greedy algorithm described in the previous section, a total of 10 iterations was required in order to meet the tolerance for the relative residual of tol = 10 −5 for all frequency points. Accordingly, a GMRes scheme has been applied for the solution at 10 out of 171 frequency points. These solutions correspond to the basis vectors used for the least squares approximations in the inner loop of Algorithm 1. In order to ensure sufficient accuracy of the basis vectors, a tolerance of gmres = 10 −8 was chosen for the GMRes scheme resulting in an elapsed time of 7.11 s. Moreover, 1562 matrix vector products were evaluated for the setup of the least squares problems resulting in 27.28 s. The solution of the least squares problems required a computational time of 2.84 s. In addition with the other parts of the algorithm that only marginally contribute to the computational effort, the total wall clock time added up to 37.32 s. This corresponds to an average of 0.218 s for the solution at each frequency point, and hence, a reduction by more than 60% compared to the conventional solution.
The convergence history of the greedy algorithm is given in Figure 1. The intermediate relative residuals after the first, sixth, and ninth iteration as well as the residual after convergence are shown. An initial frequency of 125 Hz has been chosen for the first iteration, which is the mid value of the considered frequency range. As expected, the solution at 125 Hz alone is not sufficient to approximate the solutions at other frequency points and hence only yields a decrease of the residual in its immediate vicinity. The residual after the sixth iteration exhibits six sharp minima (40, 50, 60, 125, 188, and 206 Hz), each corresponding to a GMRes solution in the outer loop and hence falling below the defined tolerance of This example indicates the effectiveness of a greedy scheme in the context of a multifrequency analysis. Choosing the frequency at which the approximation is worst as the next sample sooner or later leads to evaluation of the frequency dependent system at (or near) its resonance frequencies. In the present example, convergence was achieved rather quickly once the responses close to the resonances were added to the reduced basis. This may lead to the intuitive assumption that including responses exactly at the eigenfrequencies of the system is essential for fast convergence of the algorithm. However, excluding the frequencies 50, 100, 150, and 200 Hz leads to the choice of neighboring frequency points (eg, 49 Hz instead of 50 Hz), and the algorithm requires the exact same number of iterations for this problem. Moreover, choosing the frequency interval in a mid-frequency range, and thereby omitting the responses around the first couple of eigenfrequencies of the system, does not deteriorate the convergence of the method. For example, the solution of this problem in the frequency range from 240 to 410 Hz with frequency steps of Δf = 1 Hz required 12 iterations in order to meet the tolerance of tol = 10 −5 . The larger number of iterations compared to the low-frequency analysis (12 instead of 10) could be associated with the reduced accuracy of the solutions due to larger discretization errors. Further mathematical study is required in this regard.
Our preliminary numerical studies also indicate that the convergence behavior is rather insensitive to the choice of the initial frequency (1) . Furthermore, refining the resolution of the frequency range of interest in this example does not lead to an increase of the total number of iterations. In general, additional frequency points in a rather smooth region only affect the number of least squares solutions and not the number of actual solutions of the system. This is also reflected in Figure 2, where the solution time of the duct problem is shown for different frequency resolutions. The time required for the actual solutions of the system is unchanged while the time required for the setup and solution of the least squares problems increases linearly with the number of frequency points m; however, the rate of increase is smaller than in the conventional solution strategy. Therefore, using the greedy algorithm, the common practice among engineers of oversampling the frequency range has less impact on the overall computational effort.

Point-excited spherical shell in water
As a second example, we consider a structural acoustic interaction problem. A spherical shell made of steel is locally excited by a point force of F = 1 N and sound is radiated into the surrounding water. The geometrical parameters of the sphere and the material properties of steel and water are given in Table 1. Eight-noded quadrilateral shell finite elements based on the Reissner-Mindlin theory are employed to model the structural subdomain. Discontinuous quadrilateral boundary elements are used for the discretization of the acoustic subdomain. These boundary elements are  superparametric, ie, they are characterized by a four-noded bilinear sound pressure interpolation and a nine-noded biquadratic geometry approximation. Conforming meshes with 384 finite and boundary elements are employed corresponding to eight elements on a ∕2 arc. After forming the Schur complement (5), the system has 6924 degrees of freedom. The problem is studied at m = 138 nonuniformly spaced frequency points up to 100 Hz. A treatment for irregular frequencies is not required in the considered frequency range since the first eigenfrequency of the corresponding interior acoustic Dirichlet problem is approximately 148 Hz. 47 The absolute displacement solution at a point located at an angle of with respect to the point of excitation is shown in Figure 3 along with the analytical series solution. 48 Five resonances occur in the considered frequency range.
Due to the ill-conditioning of the Schur matrix (5), ie, condA( i ) = 10 5 · · ·10 8 , the use of an iterative solver would require preconditioning. However, due to the lack of a suitable preconditioner, and given the relatively small number of degrees of freedom, a direct solver has been employed for successively solving m = 138 linear systems requiring a wall clock time of 693.14 s. This corresponds to an average of 5.02 s for the solution at each frequency point.
The greedy algorithm required a total of 20 iterations in order to meet the tolerance for the relative residual of tol = 10 −5 for all frequency points. Consequently, a direct solver was employed at 20 out of 138 frequencies resulting in an elapsed time of 98.19 s. Moreover, 2241 matrix vector products were evaluated for the setup of the least squares problems resulting in 93.63 s. The solution of the least squares problems required a computational time of 9.36 s. In total, the wall clock time added up to 201.33 s corresponding to an average of 1.46 s for the solution at each frequency point. However, this number is not comparable to the average time that was required by the conventional strategy due to the difference in achieved accuracies.
The convergence history of the greedy algorithm is shown in Figure 4. Again, similar to the duct problem, distinct maxima in the residual emerge at the five resonance frequencies after iteration j = 12. The predefined tolerance of tol = 10 −5 is met after iteration j = 20.
Furthermore, in order to assess the quality of the least squares approximations, the relative difference with respect to the conventional solution is studied in what follows. It can be expressed by  wherex( i ) denotes the solution obtained using a direct solver and x( i ) is the solution obtained by the greedy algorithm. Trivially, = 0 at the frequencies that are chosen by the greedy algorithm for the calculation of the basis vectors. The relative differences are shown in Figure 5 for different tolerances tol of the relative residual. The respective numbers of iterations and wall clock times of the greedy algorithm are provided in Table 2. For instance, a predefined tolerance of tol = 10 −5 yields a relative difference ( i ) ≤ 10 −4 in the whole frequency range 1 , … , m . However, for tolerances smaller than tol = 10 −8 in this problem, the proposed scheme is less efficient than a conventional strategy using a direct solver for each frequency point.
These results lead to the conclusion that an a priori estimation of the condition numbers in the frequency range could provide indication of suitable tolerances tol and consequently facilitate the decision whether to use the greedy algorithm or a conventional strategy. Relatively inexpensive algorithms for estimating 1-norm condition numbers exist 49,50 and the thereby gained knowledge could also be incorporated into the choice of the next sample (10). This is certainly an issue requiring future research.

Honeycomb sandwich panel in air
Beside the case of structures submerged in water, structural acoustic interaction can also have a significant influence on the vibration response of lightweight structures in air, particularly for sandwich panels with large radiating surfaces. Clarkson and Brown were among the first to experimentally quantify the energy dissipation by virtue of sound radiation for a sandwich panel with honeycomb core. 51 They have shown that this phenomenon-commonly denoted as acoustic radiation damping-can exceed the extent of material-inherent damping by an order of magnitude. In the following, we apply the greedy algorithm for a vibro-acoustic analysis of a similar sandwich panel as the one tested by Clarkson and Brown. 51 However, we point out that the focus of the present analysis is set on the numerical analysis rather than on the acoustic radiation damping of sandwich structures.
The geometry of the six-sided panel with a central cut-out is shown in Figure 6. It is composed of an aluminum honeycomb core of 29 mm thickness and two 0.28 mm thick aluminum face sheets. The isotropic material properties of aluminum as well as the (assumed) equivalent orthotropic properties of the core are given in Table 3. In this example, additional structural damping is neglected, and hence, sound radiation is the only source of energy dissipation. The panel   Eight-noded quadrilateral shell finite elements and 20-noded hexahedral solid finite elements are employed for the discretization of the core and face sheets, respectively, resulting in a total of 22 131 displacement degrees of freedom. To model the acoustic subdomain, quadrilateral boundary elements with linear discontinuous sound pressure interpolation and a total of 1504 pressure degrees of freedom are employed. The nonconforming structural and acoustic meshes are coupled on the face sheets of the panel. In order to verify the accuracy of the discretization, the radiated sound power has been compared to a model with 10 times more degrees of freedom, resulting in a relative difference of less than 2.6% (0.11 dB) in the considered frequency range. A treatment for irregular frequencies is not required since the first eigenfrequency of the corresponding interior acoustic Dirichlet problem is approximately 6000 Hz.
The coupled system of Equations (4) is solved using a GMRes scheme with an incomplete LU factorization (ILU) for preconditioning the structural system K − 2 M. Following the conventional frequencywise strategy, a total of 8623 iterations were required for solving all the m = 196 systems individually subject to a relative tolerance of tol = 10 −5 . The total solution time including setup of the preconditioners added up to 105 minutes. The radiated sound power and the total potential energy of the vibrating panel are shown in Figure 7. While the first couple of resonances exhibit sharp peaks, the effect of acoustic radiation damping is clearly noticeable in the higher frequency range.
The greedy algorithm required a total of 43 outer iterations for the solution in the whole frequency range. The responses at the sampling frequencies were computed using the same preconditioned GMRes scheme as above with gmres = 10 −8 accounting for a wall clock time of 26 minutes. The setup and the solution of the least squares problems resulted in an elapsed time of 7 minutes. In total, the greedy algorithm required 33 minutes for the solution at all frequency points with a relative residual of less then tol = 10 −5 . This corresponds to a speedup of more than three times compared to the conventional solution.

CONCLUSION AND FUTURE WORK
A greedy reduced basis scheme has been proposed for the frequency range solution of vibro-acoustic problems. It is based on iteratively expanding the reduced basis by adding the frequency response that is currently worst approximated. The method has been applied to an interior acoustic problem as well as to coupled structural exterior acoustic problems. In all cases, convergence was reached relatively fast, and consequently, the actual solution of the system was only required at a few frequency points. The solutions at the other intermediate frequencies could be accurately approximated by linear combination using a least squares solver. Comparisons to conventional frequencywise strategies indicate the efficiency of the proposed scheme.
Although in all examples, frequencies near resonances were chosen for calculating the basis vectors during the iterations, preliminary studies have shown that explicit exclusion of eigenfrequencies from the frequency range of interest does not deteriorate convergence. Moreover, the algorithm seems to be insensitive to the choice of the initial frequency. Furthermore, refining the frequency resolution in a smooth region leads only to a marginal (if at all) increase in the total number of iterations. In general, refining the frequency resolution over a certain threshold only increases the effort required for the setup and solutions of the least squares problems. Consequently, the common practice among engineers of oversampling the frequency range has less impact on the overall computational effort than when using a conventional frequencywise strategy.
The numerical examples were carried out using the same excitation for all frequency points. While the method allows for a straightforward application of frequency-dependent right-hand sides, the success of the greedy algorithm will in these cases depend on the regularity of the response with respect to the frequency. When the responses at different frequencies are not or hardly related to each other, conventional frequencywise strategies could be more efficient. Furthermore, when considering multiple forcing vectors, adding the responses for all excitations to the reduced basis at each frequency sample could diminish the efficiency of the method. As a remedy, a truncated singular value decomposition could be applied in order to reduce the number of algebraically independent right-hand sides. Detailed studies on the use of MOR in the context of fast frequency sweeps with FEM and many forcing vectors are available in the literature. 52,53 Similar studies are also planned for the herein presented greedy reduced basis scheme.
Moreover, in a future work, a priori estimation of the condition numbers could provide indication of suitable tolerances along with an a priori error estimation and also improve the choice of basis vectors. Furthermore, the applicability of the scheme will be extended to large-scale problems by incorporating data sparse formats 26 as well as by parallelizing the algorithm. The latter could be done by solving the least squares problems in a distributed memory environment and solving the actual linear systems on a shared memory environment. Systematic benchmarking will be conducted to verify the applicability of the proposed scheme to large-scale multifrequency problems.