Model order reduction with Krylov subspaces of exterior acoustic BEM
systems
Dionysios PANAGIOTOPOULOS(1),(2)∗ , Elke DECKERS(1),(2) , Wim DESMET(1),(2)(1) Mechanical engineering department, KU Leuven(2) DMMS lab, Flanders Make

Abstract
Numerical modelling with the Boundary Element Method is often hindered due to the excessive computational
cost at higher frequencies. In that context, a suitable model order reduction scheme is proposed to alleviate
the computational cost needed for a BE fast frequency sweep acoustic analysis. First, the method utilizes akernel series expansion to remove the frequency dependency from the BEM system matrices. Subsequently, aGalerkin projection is deployed on the frequency independent matrices to reduce the size of the system. The
basis employed in the Galerkin projection is assembled by the Krylov subspaces obtained through an Arnoldi
procedure of the BEM system on a predefined frequency grid and is augmented to enable the CHIEF method
extension. The use of Krylov subspaces facilitates the definition of an error estimator that indicates the levelof the error expected by the projection of the system. The efficiency of the model order reduction scheme isvalidated on a test case.
Keywords: Model Order Reduction, BEM, Krylov subspaces recycling, Arnoldi, error estimator

1

INTRODUCTION

The Boundary Element Method [3] constitutes a widely used numerical tool in acoustic analyses. Nevertheless,the computational cost accompanying the deployment of the method renders often frequency sweep analyses
highly inefficient. The method can be partitioned into two distinct operations, namely i) the assembly of thesystem matrix and ii) the solution of the linear system of equations. In contrast with the Finite Element Method,
the system matrix is fully populated, frequency dependent and in general, it does not feature any properties suchas positive definiteness [10]. Consequently, regarding a frequency sweep analysis, the BEM necessitates assembling and solving a system of equations for each frequency for which acoustic response is desired separately.
In that context, several techniques to speed-up the calculation of acoustic response with BEM have been developed. The most established arise from the combination of the Fast Multiple Method [4] and the H-matrix[5] with BEM. Although these techniques offer substantial acceleration, they still necessitate the assembly andsolution of the system for each frequency under consideration. In that respect, several techniques have emerged
as well that aim at simplifying the multi-frequency BEM calculations by avoiding to assemble the system foreach frequency. These techniques employ approximation methods such as the discrete form interpolation method[7] or Taylor expansion for the approximation of the Green’s function [16]. Nevertheless, the effect of thesetechniques is on the one hand only limited to the assembly procedure, which is often the less computational
demanding part, and on the other hand they are often accompanied by a higher memory consumption which canbe proven critical in BEM calculations. A technique that alleviates the computational cost associated with bothprocedures of BEM was developed by Lefteriu et al. [8], where a frequency interpolation approach is combined with a dedicated Well-Conditioned Asymptotic Waveform Evaluation [15] technique for BEM. However,
the induced accelerations are only deployed sequentially, and thus the system matrix is first assembled and thenprojected to a moment-matching basis, the calculation of which might also be a cumbersome procedure.
∗ Dionysios.Panagiotopoulos@kuleuven.be7439

In this work, a novel technique for model order reduction in Boundary Element frequency sweep simulations
is introduced. The proposed technique is based upon the off-line Galerkin projection of the BEM system on abasis created by an appropriate combination of Proper Orthogonal Decomposition [9] and Krylov subspaces [2].
Specifically, the employed projection basis construction is inspired through the Krylov subspaces recycling [12],in which the Krylov subspaces produced for the jth system of a set of slowly varying systems are utilised aspart of the projection subspaces of the ( j + 1)st system. In detail, the subspaces that describe the approximatedeigenvectors (i.e. Ritz vectors) of the jth system attempt to describe as well the Ritz vectors of the ( j + 1)st
system. Analogously, in this work the projection basis is constructed by accumulating and orthogonalizing the
Krylov subspaces of BEM systems on a predefined frequency grid. CHIEF overdetermination [14] of the systemis enabled by appropriately augmenting the basis.
Furthermore, the projection is conducted in an off-line fashion. After the basis is acquired, a Taylor expansionof the frequency coupled Green’s function kernel is leveraged to approximate the BEM system as a series offrequency independent matrices. These matrices are assembled and projected gradually with respect to eachpower of the series, in order to ensure that the memory requirements do not exceed those of a single BEM
system. After the projection of all matrices of the series, the system at all frequencies is described by a seriesof frequency independent matrices that scale with the order of the projection basis. Employing a second moreelaborated reduced model, an error estimator for the first model is also defined.
The work in this paper constitutes an extension of the method presented in [11] and it is structured as follows.
In section 2 the BEM system of equations is derived and formulated as a series of frequency-decoupled matrices
by an appropriate Taylor expansion. Next, in section 3 a Galerkin projection is deployed in the context of
the Series Expansion BEM, the projection basis construction procedure is elaborated and an error estimator isdefined. In section 4 the computational gain that is offered for a benchmark model of acoustic simulations ispresented. Finally, the paper is concluded in section 5.

2

BOUNDARY ELEMENT METHOD FORMULATION WITH A SERIES EXPANSIONIn this section a direct boundary element formulation is derived through a collocational approach. Subsequently,
employing a Taylor expansion, a Series Expansion BEM (SEBEM) is presented.
2.1 BEM formulation
The representation formula (1) upon which BEM is based can be derived from the well-known homogeneous
Helmholtz equation and is given as follows:∂ G(x, y)p(x)dΓtot (x) = jρ0 ωG(x, y)un (x)dΓtot (x),Γtot
Γtot ∂ n(x)Zc(y)p(y) +Zy ∈ Ωtot ,x ∈ Γtot .

(1)

It expresses the pressure at any location of the domain Ωtot in terms of the values of acoustic pressure p andfluid velocity normal to the surface un at all boundaries Γtot of a given geometry. In detail, in (1) ω = 2π f
is the angular velocity, ρ0 is the density of the fluid, x, y are position vectors and un (x) is the normal vector
pointing outside of the domain of interest (exterior, interior). The participation factor c(y) is considered as theexterior solid angle at y and determines to what extend the specific position lies in the domain of interest.
Furthermore, the integrals in the left hand-side of (1) are the single layer potentials while in the right hand-side
the integrals represent the double layer potentials. Additionally, G(x, y) represents the Green’s function relating
any two points of the domain Ωtot and is defined as:
G(x, y) =1 e− jkr(x,y),
4π r(x, y))(2)where r(x, y) is the Euclidean distance between two locations of the domain, k = cω0 is the wavenumber and c0
the speed of sound in the fluid.

7440

Discretizing the boundary Γtot by a number of elements N, the integrals are transformed to element integralsthat are simpler to calculate. Collocating position vector y on the different nodes of the geometry, a system of
equations is obtained relating the acoustic pressures p with the normal fluid velocities un on the boundary by
G, H ∈ CN×N as follows:
Hp = Gun .

(3)

Appending the boundary conditions of the problem, expression (3) is transformed to a well-determined linear
system of equations:A( f )z( f ) = b( f ),(4)

where A ∈ CN×N . As denoted in (4), the defined linear system involves a frequency-dependency and thus, itneeds to be assembled and solved for each frequency line under consideration.
2.2 Series Expansion BEM
The Series Expansion BEM (SEBEM) is based on the Taylor expansion of the Green’s function (2) to convert
the frequency coupled part e− jkr into a power series with frequency-independent coefficients. By utilizing series
expansion to approximate the Green’s function, the frequency dependency previously lying in the kernel of thesingle and double layer potentials of equation (1) is now removed as demonstrated by:Z

1
G(x, y)dr =4πZ

e− jkr
1 Mmax (k − k0 )m
dr =
∑ m!
r
4π m=0Ze− jk0 r(− jr)m dr,r(5)

for the single layer potential, where k0 is the wavenumber corresponding to the frequency at which the Taylor
expansion is performed and Mmax the order of the expansion. Thus, although a different kernel is involved foreach coefficient of the power series, these need to be calculated only once per Taylor expansion and not foreach frequency the evaluation of acoustic response is desired. Introducing (5) and the respective expression forthe double layer potential in (3) and (4), the system takes the form of: Mmax (k − k )m

∑

m=0

0

m!


 Mmax (k − k )m0
Gm un = ∑
Hm p
m!
m=0

(6)


Mmax
(k − k0 )m
Am z( f ) = ∑
bm ,
m!
m=0

(7)

and:
 Mmax (k − k )m

∑

m=0

0

m!

respectively, where G, H are transformed to polynomial matrices with coefficients Gm , Hm ∈ CN×N .
This procedure becomes cost-efficient in case the number of frequencies of desired acoustic response is greaterthan the order of the Taylor expansion employed, i.e. N fmax > Mmax , which is often the case in frequency sweep
analyses. Nevertheless, even in these cases the SEBEM might prove to be more computationally demanding
with respect to the storage cost. In fact, it demands the simultaneous storage of Mmax square matrices thatraises the memory scaling to O(Mmax × N 2 ), while the memory cost with conventional BEM scales with acomparatively moderate O(N 2 ).

3

MODEL ORDER REDUCTION IN THE BOUNDARY ELEMENT METHOD

In this section, the model order reduction of BEM system is illustrated. First, the reduction of the BEM systemthrough a Galerkin projection is elaborated. Subsequently, the projection basis employed for the model reduction
is derived and finally, an error estimator for the technique is proposed.

7441

3.1 Galerkin projection in SEBEM
Bringing the BEM system into the form of equations (6, 7), offers the benefit that the matrices Gm , Hm are
frequency independent and enables the accurate approximation of systems for the whole frequency region the
Taylor expansion is valid. Exploiting this benefit, a Galerkin projection is deployed to succeed in the order reduction of the matrix. The Galerkin projection is based on the approximation of the unknown acoustic variables
by a linear combination of linearly independent vectors, defined by the column vectors of V` ∈ CN×` :z ≈ ẑ = V` z` .

(8)

Substituting equation (8) in (4) yields a right multiplication of the system matrix A.
Left-multiplying with another projection matrix WT , with W ∈ CN×` , the system is brought into a square formof dimension ` by an oblique projection. Enforcing W = V` an orthogonal or Galerkin projection is obtained.
Employing this projection after the assembly of the BEM system coincides with the procedure that is followed
by indirect solvers, such as GMRES [13], where different projection matrices are constructed for each system.
On the contrary, employing this technique in the SEBEM framework allows for a simultaneous projection of all
the systems that are approximated by the Taylor expansion, as follows: Mmax (k − k )m

∑

m=0

0

m!


Mmax
(k − k0 )m ∗V∗` Am V` z( f ) = ∑V` b m ,
m!
m=0

(9)

and is performed after the assembly of the derivative matrices Am . Projecting each derivative matrix right after
its assembly, ensures that the memory requirements remain within the limits of the conventional BEM, as thecomputational cost in terms of memory scales more favourably with O(Mmax × `2 ).
3.2 Krylov projection basis
Deploying the Galerkin projection within the BEM framework, requires that an appropriate basis be constructed. The recently presented concept of Krylov subspace recycling [12] provides the respective backgroundto create a representative basis. Specifically, the Krylov subspaces produced for a set of master frequencies
F = [ f1 , f2 , . . . , fL ] are recycled for the approximation of the acoustic response in the whole frequency range under consideration. The Krylov subspaces of order s at a specific frequency fi can be generated by the Arnoldi
algorithm [1] and are represented as:
Ksfi (A fi , b fi ) = span{b fi , A fi b fi , A2fi , b fi , . . . , As−1
fi b fi }.

(10)

Accumulating the Krylov subspaces for the set of master frequencies, a basis that spans the union of the Krylov
subspaces of all master frequencies is constructed as:
Ktot = span{Ksf1 ∪ Ksf2 ∪ · · · ∪ KsfL }.

(11)

Practically, the projection matrix is constructed by accumulating the vectors produced by the Arnoldi algorithm
for each master frequency and subsequently, orthogonalizing and truncating through a Singular Value Decomposition procedure for the subspaces of the highest energy. The basis that is finally constructed, proves to be
representative enough as it attempts to approximate the response of a system on a given frequency fξ witha subspace doubly enriched. On the one hand, it contains the subspaces of the Ritz vectors of neighbouring
frequencies and on the other hand, subspaces of Ritz vectors of multiple frequencies e.g. fη = 2 fξ , that provide
flexibility to account for the local variations of the surface variables p and un .
3.3 CHIEF method incorporation
In exterior problems simulated with the direct BEM the effect of fictitious resonances is apparent, resulting ininstabilities in the yielded frequency responses. This effect is also apparent in combination with the proposed

7442

MOR technique as documented in [11]. To mitigate these instabilities the CHIEF method [14] can be utilized.
Nevertheless, deploying the CHIEF method, the system of equations assembled in equation 4 becomes overdetermined. Thus, the Galerkin projection as described in section 3.1 is not possible. Instead, to enable the modelorder reduction of BEM systems that are augmented with CHIEF points entries ACH ( f ) ∈ CN×q , an oblique projection is employed. The basis utilized is augmented by the same subspaces that are defined by the coefficientsof ACH at the master different frequencies as follows:
f1
f2
fL
∪ ACH
· · · ∪ ACH},
Ktot = span{Ksf1 ∪ Ksf2 ∪ · · · ∪ KsfL ∪ ACH(12)

maintaining q in low levels to ensure that the basis does not grow significantly in dimensions. The oblique(N+q)×N :projection is then conducted as follows for each derivative matrix ACH
m ∈C
ACH
m,red =



V∗aug
0

0
I



Am
ACH
m





Vaug



.

(13)

Finally, projecting all the derivative matrices, a power series with coefficients unequally dimensioned matrices(l+q)×l is obtained.
ACHm,red ∈ C
3.4 Error estimator
Krylov subspaces permit also the derivation of an error estimator. This is based on the fact that employing
additional subspaces contributes to cover more space related to the approximated eigenvectors of each system.
Namely, the new subspace succeeds in describing better and more Ritz vectors for each system.
The relative error is defined as the normalized Euclidean norm of the difference between the true solution z
and the approximated solution ẑs of the system. Considering a sufficiently higher order of Krylov subspaces,
enables the exploitation of the triangle inequality to derive an error estimator as follows:εs =|ẑs − z| |ẑs − ẑs+tn + ẑq+tn − z| |ẑs − ẑs+tn | |ress+tn |=≤+|ẑs ||ẑs ||ẑs ||ẑs |(14)

Due to the fact that Krylov subspaces of higher order s +tn contain the subspaces of lower order s, the last termof expression (14) can be assumed to be negligible and be truncated. Thus, an error estimator εˆs is obtained,
as follows:εˆs =|ẑs − ẑs+tn |.
|ẑs |(15)

The definition of the error estimator implies that two Galerkin projections are necessitated in order to assess
the quality of a reduced model. Specifically, the Taylor expansion matrices are projected on the target and ona slightly more detailed basis and the difference of the resulting surface distribution of the acoustic variables
serves as a good indication of the true error. As the computational cost of solving a reduced model scales witha function of the reduced order, the total cost induced is still considerably lower than solving the full model.

4

NUMERICAL ASSESSMENT

In this section a numerical evaluation of the proposed method is performed by deploying it to the well-known
cat’s eye benchmark case. The performance of the method is assessed both in terms of accuracy and computational efficiency.
4.1 Numerical example
To assess the proposed model reduction technique, the well-documented cat’s eye model is utilized [6], as
illustrated in Figure 1 and acoustic pressure is evaluated at locations R1 = [−1, −1, −1] and R2 = [0.1, 0.1, 0.1].

7443

un = 1m/s

un = 0m/s

Figure 1. Cat’s eye sketch

Figure 2. Sound Pressure Level at Rec. 1 and 2 yielded by BEM and MOR; top figure, Relative error and error
estimate of surface pressure; bottom figure; Cat’s eye model
The model contains 3114 nodes and equal number of degrees of freedom and its validity extends till 800Hz
for the requirement of 6 elements per wavelength. Selecting a refined grid of frequencies with increment of
1Hz, the computational requirements employing conventional BEM rise to assembling and solving 751 discrete
systems for the range of interest F = [50, 800]Hz.
On the contrary, employing the proposed technique both the number of required system assemblies is reduced
and the solution of each system is accelerated. Specifically, apart from assembling 50 system derivative matrices
selecting a Taylor expansion at 400Hz, an additional assembly of 16 system assemblies is required accounting
for the master frequencies, which are defined as Fmaster = [50, 100, . . . , 750, 800]. Thus, the total number of fullsystem assemblies rises to a total of 66 matrices.
The projection basis is assembled by collecting the Krylov subspaces of 15th order produced by the systemof each master frequency. In order to better describe the solution of the overdetermined systems of equations
that are constructed with the CHIEF method to avoid spurious resonances, the subspaces are augmented by theGreen’s function of a limited number of CHIEF points. Thus, enriching the basis with the subspaces of the
CHIEF points’ Green’s function, the resulting system of equations has dimensions 305 × 294, approximately

7444

Figure 3. Computational resources comparison measured in CPU time and memory required; Cat’s eye model
10% of the size of the original system.
In Figure 2, it is illustrated that the relative error induced due to the use of the proposed technique is moderate
and within acceptable limits. The local increase of the relative errors can be associated with the resonance frequency of the interior cavity. However, due to the addition of the CHIEF points, the respective sound pressurelevels are not heavily affected and thus, the resulting sound pressure levels demonstrate the behaviour documented in the literature. The error estimator considering tn = 4 additional subspaces for each master frequency,
observes closely the curve of the true relative error especially for the high frequencies, where increasing theorder of the subspaces offers a greater improvement with respect to the coverage to the solution space than in
the lower frequencies.
In Figure 3, a comparison of the computational resources required by BEM, SEBEM and the proposed technique(MOR) is performed. Due to the small size of the system, the solution CPU time is negligible in comparison
to the assembly time for all three methods and thus, the advantage of accelerating the respective procedure is
not pronounced in absolute terms. Due to this fact, the proposed technique necessitates approximately the same
CPU time as the SEBEM. Nevertheless in relative terms, comparing the CPU time required for the solutionof the system, the MOR technique is more competitive as it scales with one order of magnitude lower than
SEBEM and BEM, namely in seconds in comparison to minutes. Finally, in the MOR technique the memoryadvantage is also apparent, as only the storage of the reduced Taylor matrices is required.

5

CONCLUSIONSIn this work, a model order reduction technique for acoustic BEM is presented incorporating the CHIEF method.
The proposed method is based on the Taylor expansion of the Green’s function kernel as a frequency decouplingstrategy and a Galerkin projection is deployed for the order reduction of the resulting polynomial system. Thus,
the proposed technique accelerates both the assembly and the solution of the system without requiring additional
memory. The projection basis is assembled levering the Krylov subspaces created for systems on a predefined
grid of master frequencies. Exploiting the Krylov subspace properties an error estimator is defined to assess thequality of the reduced model. Finally in this paper, after being introduced, the method is analysed in terms ofcomputational cost and validated on the previous documented cat’s eye model.

ACKNOWLEDGEMENTS
The research of D. Panagiotopoulos is funded by an Early Stage Researcher grant within the European Project
PBNv2 Marie Curie Initial Training Network (GA 721615) and the research of E. Deckers by a grant from the7445

Research Foundation – Flanders (FWO). The Research Fund KU Leuven is also gratefully acknowledged for its
support.

REFERENCES[1] W. E. Arnoldi. The principle of minimized iterations in the solution of the matrix eigenvalue problem.
Quarterly of applied mathematics, 9(1):17–29, 1951.
[2] Z. Bai. Krylov subspace techniques for reduced-order modeling of large-scale dynamical systems. Applied
numerical mathematics, 43(1-2):9–44, 2002.
[3] C. Brebbia and R. Butterfield. Formal equivalence of direct and indirect boundary element methods. Applied Mathematical Modelling, 2(2):132 – 134, 1978.
[4] M. Fischer, U. Gauger, and L. Gaul. A multipole galerkin boundary element method for acoustics. Engineering analysis with boundary elements, 28(2):155–162, 2004.
[5] W. Hackbusch, B. Khoromskij, and S. A. Sauter. On h 2-matrices. In Lectures on applied mathematics,
pages 9–29. Springer, 2000.
[6] M. Hornikx, M. Kaltenbacher, and S. Marburg. A platform for benchmark cases in computational acoustics.
Acta Acustica united with Acustica, 101(4):811–820, 2015.
[7] S. M. Kirkup and D. Henwood. Methods for speeding up the boundary element solution of acoustic
radiation problems. Journal of vibration and acoustics, 114(3):374–380, 1992.
[8] S. Lefteriu, M. S. Lenzi, H. Beriot, M. Tournour, and W. Desmet. Fast frequency sweep method for indirect boundary element models arising in acoustics. Engineering Analysis with Boundary Elements, 69:32–45,
2016.
[9] Y. Liang, W. Lin, H. Lee, S. Lim, K. Lee, and H. Sun. Proper orthogonal decomposition and its
applications–part ii: Model reduction for mems dynamical analysis. Journal of Sound and Vibration,
256(3):515–532, 2002.
[10] S. Marburg and B. Nolte. Computational acoustics of noise propagation in fluids: finite and boundary element
methods, volume 578. Springer, 2008.
[11] D. Panagiotopoulos, E. Deckers, and W. Desmet. Krylov subspaces recycling based model order reduction
for acoustic bem systems and an error estimator. Computer Methods in Applied Mechanics and Engineering,(submitted), 2019.
[12] M. L. Parks, E. De Sturler, G. Mackey, D. D. Johnson, and S. Maiti. Recycling krylov subspaces forsequences of linear systems. SIAM Journal on Scientific Computing, 28(5):1651–1674, 2006.
[13] Y. Saad and M. H. Schultz. Gmres: A generalized minimal residual algorithm for solving nonsymmetriclinear systems. SIAM Journal on scientific and statistical computing, 7(3):856–869, 1986.
[14] H. A. Schenck. Improved integral formulation for acoustic radiation problems. The journal of the acousticalsociety of America, 44(1):41–58, 1968.
[15] R. D. Slone, R. Lee, and J.-F. Lee. Broadband model order reduction of polynomial matrix equations
using single-point well-conditioned asymptotic waveform evaluation: derivations and theory. International
Journal for Numerical Methods in Engineering, 58(15):2325–2342, 2003.
[16] Z. Wang, Z. Zhao, Z. Liu, and Q. Huang. A method for multi-frequency calculation of boundary integral
equation in acoustics based on series expansion. Applied Acoustics, 70(3):459–468, 2009.

7446