Accelerated sound propagation simulations using an error-free Fourier method coupled with the spectral element method

Nikolas Borrel-Jensen 1

Acoustic Technology, Department of Electrical and Photonics Engineering, Technical University of Denmark, 2800 Kongens Lyngby, Denmark

Allan P. Engsig-Karup 2

Department of Applied Mathematics and Computer Science, Technical University of Denmark, 2800 Kongens Lyngby, Denmark

Maarten Hornikx 3

Department of the Built Environment, Eindhoven University of Technology, 5612 AZ Eindhoven, The Netherlands

Cheol-Ho Jeong 4

Acoustic Technology, Department of Electrical and Photonics Engineering, Technical University of Denmark, 2800 Kongens Lyngby, Denmark

ABSTRACT Simulating acoustics e ffi ciently and accurately using numerical methods has been an active research area for the last decades and has applications in computer games, VR / AR, and architectural design. However, their extensive computation time makes these methods challenging for large scenes and broad frequency ranges. This work attempts to accelerate the simulations using rectangular decomposition, enabling error-free propagation in the bulk of the domain consisting of air. We exploit the analytical solution to the wave equation in rectangular domains calculated using the Fast Fourier Transform with near-optimal spatial discretization satisfying the Nyquist criterium. Coupling with the spectral element method near the boundaries results in a method capable of handling complex geometries with realistic boundaries, though with the caveat that additional errors and computational overhead may result from the interface. This paper will investigate the accuracy and e ffi ciency of the proposed domain decomposition method compared to a spectral element implementation running in the entire domain and the results in 1D indicate an 18 times speedup factor for relative errors below 9%.

1 nibor@dtu.dk

2 apek@dtu.dk

3 m.c.j.hornikx@tue.nl

4 chje@dtu.dk

a slaty. inter.noise 21-24 AUGUST SCOTTISH EVENT CAMPUS O ¥, ? GLASGOW

1. INTRODUCTION

Simulating acoustics in virtual spaces is an active research topic, with applications in computer games, VR / AR, and building acoustics. Using numerical methods to solve the underlying partial di ff erential equations naturally takes the wave phenomena into account. However, it comes with the cost of being computational demanding, especially when simulating broad frequency ranges and large domains. To overcome the demanding computations, we will develop an e ffi cient and accurate domain decomposition method coupling an e ffi cient and error-free Fourier method in rectangular domains with an spectral element method applied at the boundaries handling complex geometries and impedance boundaries. The most used numerical methods for simulating acoustics are the finite element methods (FEM) [12], spectral element methods (SEM) [14,19], DG-FEM [11], finite-di ff erence time-domain methods (FDTD) [5], boundary element methods [6], the Wave-Based Method (WBM) [17], and Pseudo-spectral Fourier methods [8]. Domain decomposition methods (DDM) [3] are widely used approaches where the domain is split into many partitions. We will use DDMs to divide the domain into simpler partitions where specialized methods can be applied independently depending on the properties of the domain. Our approach is inspired by Raghuvanshi et al. [10,16], where an adaptive rectangular decomposition (ARD) method was proposed, exploiting the well-known analytical solution to the wave-equation in rectangular domains coupled with the FDTD method near the boundaries. The ARD method has error-free propagation within the rectangular domains and can be computed very e ffi ciently using the Fast Fourier Transform (FFT). Only numerical errors are introduced at the interfaces. We follow the same approach, but instead of coupling with FDTD, we couple the Fourier method (FM) with SEM, still using FDTD for interface handling though. The SEM is advantageous over FDTD as it handles complex geometries and impedance boundaries [14] in a robust manner. By coupling the e ffi cient Fourier method running in rectangular domains, and the SEM near the boundaries, we can speed up the overall computation time for large scenes, where the ratio between air and boundaries is large. Related domain decomposition approaches have been taken to achieve similar goals. In [13], the time-domain Pseudo-Spectral element method was coupled with the DG-FEM method near the boundaries. In [18] the Wave-Based method is applied in convex domains and coupled with the second-order SEM in domains requiring geometrical flexibility. We propose a method consisting of three parts: 1) an error-free Fourier method (FM) running in rectangular domains, 2) the SEM for modeling complex geometries and realistic impedance boundaries, and 3) an FDTD scheme for interface handling. In Figure 1a the domain decomposition is illustrated.

Figure 1: Overview of the domain decomposition methods coupling the Fourier method in rectangular partitions with the spectral element method near the boundaries.

SEM

interface

FOURIER

METHOD

interface

interface

interface

interface

FOURIER METHOD

SEM

interface

interface

(b) Time step from n to n + 1 with Neumann boundary conditions at the interface for each of the methods running in the separate partitions. Note, that the methods in the partitions are running independently and compensation of the reflections at the interface is done after each time step.

(a) A domain decomposed into rectangular partitions running the error-free FM coupled with the SEM near the boundaries.

time step n> n+1 Neumann p(n+t) interface handling — pint) ---> after interface handling

2. GOVERNING EQUATIONS

The wave equation is

∂ 2 p ( x , t )

∂ t 2 − c 2 ∇ p ( x , t ) = F ( x , t ) , t ∈ R + , x ∈ R N . (1)

where p is the pressure (Pa), t is the time (s) and c is the speed of sound in air (m / s) and F ( x , t ) is a forcing function. An initial conditions can be used instead of a forcing function satisfied by using e.g. a Gaussian source for the pressure part and setting the velocity equal to zero

! 2  , ∂ p ( x , t = 0)

p ( x , t = 0) = exp  − x − x 0

∂ t = 0 , x 0 ∈ R N , (2)

σ 0

with σ 0 being the width of the pulse determining the frequencies to span and boundary conditions

∂ p ( x , t )

∂ n = v , x ∈ Γ v , (3)

where v is the enforced velocity at the boundary Γ v , and n is the normal pointing outwards from the boundary Γ . Equation (3) is the Neumann boundary conditions and is applied at the interface for coupling domains – as explained in Section 4 – with v = 0.

3. THE FOURIER METHOD

3.1. The analytical solution It can be shown by using separation of variables [1] (p. 155-157), that any 1D pressure field p ( x , t ) in rectangular domains with Neumann boundary conditions can be represented in the form of a general series representation

N − 1 X

i = 0 M i ( t ) Φ i ( x ) , where M i ( t ) = ae jck i t + be − jck i t , k i = π i

l , (4)

p ( x , t ) =

and i is the mode, N is the maximum mode to include depen din g on the required frequency range, M i ( t ) is the time-varying mode coe ffi cients [9, 15], j = √

− 1 is the complex number, k i is the wavenumber, and Φ i ( x ) = cos ( k i x ) are the eigenfunctions of the Laplacian for a rectangular domain. The time constants a and b depend on the initial conditions.

3.2. The discrete formulation We can interpret the analytical solution (4) on a discrete uniform grid x i with the highest wavenumber being spatially sampled at the Nyquist rate. Assuming that the signal is properly band- limited and that enough modes N are included to capture the band-limited signal, the discretization introduces no numerical errors. In the discrete interpretation, Equation (4) is similar to the inverse Discrete Cosine Transform, with ϕ i being the cosine basis vectors. Hence, we can convert from modal coe ffi cients M to pressure values p as p ( t ) = iDCT( M ( t )) , (5)

where p is the N × 1 pressure vector and M ( t ) is the N × 1 vector including all modes for time t . Reinterpreting the wave equation ∂ 2 p

∂ t 2 − c 2 ∂ 2 p

∂ x 2 = F in a discrete setting for the spatial dimensions and taking the cosine transform of both sides of the equality sign yields

∂ 2

∂ t 2 M ( t ) + c 2 k 2 ⊙ M ( t ) = DCT( F ( t )) , (6)

where F ( t ) is a N × 1 vector with pressure values of the forcing function in all spatial nodes at time t , M ( t ) is a N × 1 vector with mode coe ffi cients corresponding to time t , k is a N × 1 vector including the wave numbers, and ⊙ is the Hadamard product operator. Disregarding the forcing term, the above equation describes a set of independent simple harmonic oscillators, each vibrating with its characteristic frequency ω i = ck i . In the case where the forcing term is a constant transformed into mode-space ˜ F , we can simply solve the equation ∂ 2

∂ t 2 M i + ω 2 i M i − ˜ F i = 0, but since the forcing term changes with time, we need to derive a temporal update scheme for this. We will assume that the forcing term is constant over a time-step ∆ t , which is satisfied when proper Nyquist sampling is applied. The forcing term may be transformed into mode-space as ˜ F ( n ∆ t ) ≜ DCT( F ( n ∆ t )). We will use the second-order centered di ff erence multiplied by the term ω 2 ∆ t 2 2(1 − cos( ω ∆ t )) originating from the solution to the simple harmonic oscillator, obtaining the update scheme [16]

M n + 1 i = 2 M n i cos( ω i ∆ t ) − M n − 1 i + 2 ˜ F n

ω 2 i (1 − cos( ω i ∆ t )) . (7)

3.3. Discrete Fourier Transform for imposing Neumann boundary conditions The formulation in the analytical formulation in (4) is real and hence, to ensure a real-valued frequency spectrum after applying the Discrete Fourier Transform, an even Fourier continuation p k = p − k and periodicity 5 are required. Mirroring the samples around a boundary satifies the zero derivative Neumann boundary condition in that point. The extended time signal ˆ p k of lenght (2 N − 2) is then

 p k , for 0 ≤ k < N ,

ˆ p k

(8)

p 2 N − k , for N ≤ k < 2 N − 1 .

The single-sided frequency spectrum M i consisting of N bins is used for time-stepping in Equation (7). Again, but now in the Fourier domain, the even and periodic double-sided spectrum ˆ M i of size (2 N − 2) is reconstructed and the inverse Fourier Transform is applied to retrieve the time pressures back as

 M i , for 0 ≤ i < N ,

ˆ M i

(9)

M 2 N − i , for N ≤ i ≤ 2 N − 1 .

Extracting the first N pressures values from F − 1 ( ˆ M i ) gives us the final result. Hence, the DCT and iDCT used in Section 3 and the remaining sections are referring to the definitions below

DCT  p k  = Re ( F [ ˆ p k ]) ,

iDCT[ M i ] = Re  F − 1 [ ˆ M i ]  . (10)

4. INTERFACE HANDLING

We will derive the interface scheme for communicating pressure values between partitions following the same approach as Raghuvanshi [15,16] to couple the FM in rectangular partitions with the SEM near the boundaries. For clarity of presentation, we will consider the 1D wave equation using the second-order centered di ff erences. The wave equation (1) discretized in space can be written as ∂ 2 P

∂ t 2 − KP = F ( t ), where K represents the Laplacian operator at each node as a Discrete Laplacian Matrix and F is the forcing term at each node. As a result of spatial discretization, K is transformed

5 if we have a signal ‘abcd’, then the even extension would be ‘abcdcb’.

into a symmetric matrix

...





1 − 2 1

K = c 2

1 − 2 1

. (11)

∆ x 2

1 − 2 1

1 − 2 1 ...

Assume that the grid consists of N nodes. Care must be taken at the boundary, since calculating the pressure value p N at the right-most node N would require the neighboring pressure value p N + 1 located outside the grid for a Neumann boundary condition imposed by mirroring the pressure half- way between the nodes or at the boundary node. We will denote the pressure values at these ghost nodes as the residual . Eliminating the ghost nodes happens by discretizing the BCs and inserting these expressions into the scheme. For example, for a right boundary node x N , the residual is added to the inner partition as

Neumann N + 1 / 2 p ( x N − 1 ) − p ( x N )

p ′′ N + 1 / 2 ( x N ) = p ( x N − 1 ) − 2 p ( x N ) + p ( x N + 1 )

∆ x 2 ,

∆ x 2

(12)

Neumann N = 2 p ( x N − 1 ) − 2 p ( x N )

p ′′ N ( x N ) = p ( x N − 1 ) − 2 p ( x N ) + p ( x N + 1 )

∆ x 2 .

∆ x 2

These observations can be used to formulate the interface handling between two separate partitions. Suppose we wish to split the domain into two equal partitions with N nodes each, such that each might be treated independently. P is a vector of length 2 N . This can be done by decoupling K into a block diagonal form while accounting properly for the o ff -diagonal entries K = A + C , where the decoupled matrix A and the coupling matrix C are given by (interface located at N + 1 / 2)

...









0

...

1 − 2 1

1 − 1

A = c 2

, C = c 2

1 0

. (13)

∆ x 2

∆ x 2

− 1 1

0 1

...

1 − 2 1 ...

0

The pressure is first updated using A , and the residual part not taken into account is then computed using C . The coupling matrix C can intuitively be seen as a communication step to pass pressures between two partitions, and the action of C is to compute the gradient of the pressure at the interface (scaled by some factor). Rewriting in terms of the decoupled and coupled matrices yields ∂ 2 P

∂ t 2 − AP = F ( t ) + CP , where the coupling term CP can be seen as an additional source term accounting for the alignment of pressures between partitions as ˆ F ( t ) = F ( t ) + CP . We can formulate a unified framework for communicating pressures between partitions running any method. The boundary condition will be imposed at the pressure node N (not N + 1 / 2 as above) since this will ease the implementation for the SEM where boundaries are explicitly defined at the nodal points. Without lack of generality, assume that two independent (2,2) FDTD update schemes are running in each partition and denote the pressures in partition 1 as p 1 , i and partition 2 as p 2 , i with subscripts denoting the partition and corresponding node index, respectively. Then, the interface communication can be handled as follows:

Figure 2: Interface handling with twice refined spatial and temporal resolutions in the right partition compared to the left partition. (a) Neumann reflections present at time n + 1 and n + 1 / 2 with spatial interpolation in partition 1. (b) Interface handling for time n + 1 and n + 1 / 2 by using the corresponding ∆ x 1 , ∆ t 1 , ∆ x 2 , ∆ t 2 values for the FDTD scheme in Equation (14)-(16). (c) Temporal interpolation in the coarse domain for calculating n + 1. (d) Interface handling for the fine domain (only) at time n + 1.

n+ 1

n+1

n-1/2 n n+1/2

n-1/2 n n+1/2

n

n

n- 1

n- 1

(a)

(b)

n+1 n+1

n+1

n+ 1

n+1/2

n+1/2

n n+1/2

n

n

n-1

n-1

Time interpolation

(c)

(d)

1. Calculate the pressures in each domain as completely independent partitions with Neumann boundary condition at the interface

p ( n + 1) 1 , N ← c 2 ∆ t 2 1 ∆ x 2 1

 2 p ( n ) 1 , N − 1 − 2 p ( n ) 1 , N  + 2 p ( n ) 1 , N − p ( n − 1) 1 , N + F ( n ) 1 , N ,

 − 2 p ( n ) 2 , 1 + 2 p ( n ) 2 , 2  + 2 p ( n ) 2 , 1 − p ( n − 1) 2 , 1 + F ( n ) 2 , 1 . (14)

p ( n + 1) 2 , 1 ← c 2 ∆ t 2 2 ∆ x 2 2

2. Remove the residual part at time n corresponding to the reflected pressures

p ( n + 1) 1 , N ← p ( n + 1) 1 , N − c 2 ∆ t 2 1 ∆ x 2 1

 p ( n ) 1 , N − 1  , p ( n + 1) 2 , 1 ← p ( n + 1) 2 , 1 − c 2 ∆ t 2 2 ∆ x 2 2

 p ( n ) 2 , 2  . (15)

3. Transfer the removed residual part at time n to the neighboring domains

p ( n + 1) 1 , N ← p ( n + 1) 1 , N + c 2 ∆ t 2 1 ∆ x 2 1 ·  p ( n ) 2 , 2  , p ( n + 1) 2 , 1 ← p ( n + 1) 2 , 1 + c 2 ∆ t 2 2 ∆ x 2 2

 p ( n ) 1 , N − 1  . (16)

Separating the pressure update scheme inside the cavity from the interface handling scheme makes it clear that we can transfer pressures between partitions running any schemes. The sixth-order update scheme used in this work can be derived in a similar manner. The spatial and temporal resolutions will di ff er between partitions if the e ffi ciency of the Fourier method is to be exploited and therefore piece-wise cubic Hermite interpolation is used in both space and time [2]. The procedure is depicted in Figure 2.

Figure 3: Domain decomposition coupling FM with SEM by introducing a first-order SEM layer to stabilize the scheme.

0 L0

L1

L

x

P=1 Fourier

SEM

P>1 SEM

L0 L0

0

L1

L1

L

SEM coupling FD coupling

ε

4.1. Interfacing with the spectral element method Since the interface handling is independent of the methods running in the partitions, coupling the FM with the SEM can be implemented using the framework described. However, due to the non-physical behavior caused by the Neumann condition at the interface at each time-step before compensating, non-smooth second derivatives are introduced. This type of behavior is known as ‘shocks’ in the literature, and it is well-known to introduce challenges for the SEM. An investigation has been made by comparing the Laplacian term from Equation (1) for the SEM and FDTD methods

SEM Laplacian

(SEM) p ( n + 1) = 2 p ( n ) − p ( n − 1)

− c 2 ∆ t 2 M − 1  S p ( n ) 

+ ∆ t 2 f ( n ) , (17)

FDTD Laplacian

+ c 2 ∆ t 2

(FDTD) p ( n + 1) = 2 p ( n ) − p ( n − 1)

∆ x 2 K p ( n ) + ∆ t 2 f ( n ) , (18)

where M is the mass matrix and S is the sti ff ness matrix [4,7]. We have noted that noticable bigger errors are introduced for higher-order SEM compared to first-order SEM. A simple remedy was to add an SEM layer of first-order polynomials near the interface as illustrated in Figure 3.

5. NUMERICAL EXPERIMENTS

We will investigate the accuracy and e ffi ciency of the domain decomposition method coupling the FM with the SEM. All experiments are done in a 1D domain with f max = 1000 Hz with the di ff erentiated Gaussian pulse injected halfway into the FM partition 1. In the following, we will use the notation FM{ppw} and SEM{ppw}, i.e., FM4-SEM8 would be the FM-SEM coupling with spatial resolutions of four and eight points per wavelength (ppw) per partition, respectively. To assess the accuracy of the method, we will separately consider interface errors and errors over time at a receiver position. The interface errors are investigated by measuring the reflected pressures in the vicinity left to the interface for a left to right traveling wave at time t refl = 0 . 0135 s corresponding to the wave having made a single pass through the interface. The errors over time at the receiver position is compared against a FM reference solution for a total running time t max = 0 . 2 s using di ff erent measures.

5.1. Interface errors The domain length in this experiment is 10 m split into two partitions of 5 m each. We measure the interface pressure errors in dB as

! dB , (19)

max( | p interface | )

ϵ interface = 20 log 10

| p inc |

where p interface = { p ( t , x ) | x ∈ [3 . 5 , 5] m , t = 0 . 0135 s } are the reflected pressure values in the vicinity of the interface after a single interface traversal and has been chosen such that no other wavefronts

Table 1: Errors for the FM-SEM coupled using a sixth-order FDTD scheme for di ff erent spatial resolutions determined by the points per wavelength, keeping the frequency range fixed. (a) FM- SEM interface errors,(b) Relative error ϵ rel and maximum error ϵ ∞ calculated by comparing with an exact Fourier reference solution.

Interface errors at t = 0 . 0135 s

SEM 3 4 6 8 9 12

FM

3 -25 dB -26 dB † -28 dB -26 dB † -26 dB -25 dB

4 - -34 dB -35 dB † -36 dB -32 dB † -32 dB (a) † : ∆ t 1 = ∆ t 2 , since the spatial resolution is not a multiple of each other.

Errors at t = 0 . 2 s

FM ↔ SEM [ppw] ϵ rel ϵ ∞ 4 ↔ 4 - 0.2528

4 ↔ 8 9.1% 0.0811

4 ↔ 12 5.7% 0.0562 (b)

Figure 4: FM4-SEM8 running for t max = 0 . 2 s. From left to right: 1) Wave propagation in the full domain at t = 0 . 2 s, 2) impulse response at receiver position x = 6 . 0 m, 3) L 1 IR error, 4) transfer function at x = 6 . 0 m.

FM4-SEM8

1

0.5

FM-SEM reference receiver

FM-SEM reference

0.5

0

0

-0.5

-0.5

-1

0 2 4 6 8 10 x [m]

0 0.05 0.1 0.15 0.2 time [s]

-40

0.08

-50

0.06

-60

0.04

-70

0.02

FM-SEM reference

-80

0

-90

10 2 10 3

0 0.05 0.1 0.15 0.2 time [s]

Frequency [Hz]

are overlapping. The pressure values before traversing the interface are normalized by the incident pressure p inc in the domain. ϵ interface = −∞ corresponds to zero pressures (no interface error) being reflected. The interface errors for the FM-SEM coupling are summarized in Table 1a. The errors are measured for combinations of spatial resolutions in the partitions denoted by the number of ppw. Overall, we see that using three ppw for the FM shows large errors indicating that four ppw is the lowest resolution possible for reasonable accuracies. The FM can handle two ppw corresponding to the Nyquist criterium but is limited by the sixth-order FD scheme at the interface. FM4-SEM8 gives the lowest interface error of -36 dB.

5.2. Overall accuracy Another investigation of the accuracy is to compare the simulation with a reference solution. The error of the impulse responses (IRs) at receiver position x = 6 . 0 m is calculated in the relative error ϵ rel and the maximum error ϵ ∞ norms

T − 1 X

| p sim ( t n , x ) − p ref ( t n , x ) |

ϵ rel ( x ) = 1

| p ref ( t n , x ) | , ϵ ∞ ( x ) = max n = 0 , 1 ,..., T − 1 | p sim ( t n , x ) − p ref ( t n , x ) | . (20)

T

n = 0

The simulated pressures p sim ( t n , x ) and the reference pressures p ref ( t n , x ) correspond to the IRs at the receiver position x for the discrete time steps t n , and T = ⌈ t max / ∆ t ⌉ is the total number of time steps. The results are summarized in Table 1b. First, we notice that FM4-SEM4 gives the largest errors ϵ ∞ = 0 . 25 and very big relative errors (not included in the table). In contrary, coupling with SEM8 or SEM12 drastically improves the accuracy with errors of ϵ ∞ < 0 . 08 and ϵ rel < 9%. These results are not all consistent with the interface errors; notice for example, that the FM4-SEM8 coupling (-36 dB) compared to FM4-SEM12 (-32 dB) has 4 dB lower reflection errors, which could indicate that the SEM is introducing numerical dispersion errors also contributing to the overall error. Using other time integration schemes, such as the Runge-Kutta method, should decrease the numerical SEM errors. In Figure 4, the FM4-SEM8 is plotted against a reference solution showing the domain pressure at t max = 0 . 2 s, the IR and the transfer function (TF) at receiver position x = 6 . 0 m, and the corresponding L 1 errors over time for the IR. We see a good match between the simulation and the reference, though some small noticeable pressure perturbations are primarily due to interface errors.

5.3. Convergence In Figure 5 the convergence is plotted pair-wise for FM{2,4,8,16,32}-SEM{4,8,16,32,64} with 1) fixed temporal resolution, 2) individual temporal resolution satisfying the Courant-Friedrichs-Lewy condition, and 3) pair-wise for FM{16,32}-SEM{16-2,32-4} where the SEM P = 1 interface layer is oversampled eight times compared to the main SEM with fixed temporal resolution. When the P = 1 layer is finely oversampled, the forth-order spectral convergence is preserved, indicating that the interface errors are converging at the same rate.

5.4. E ffi ciency The theoretical speedup in 3D when running FM instead of SEM in the full domain is 2 × r 3 , where r is the spatial resolution factor between the FM and the SEM running in the two partitions, resulting in a 16 × speedup for FM4-SEM8 and 54 × speedup for FM4-SEM12. The factor of 2 stems from the time resolution being twice as coarse for the FM. On top of that, the SEM time complexity for solving the system of equations consisting of sparse band matrices can be done in O ( q 2 n ) + O ( qn ) in time using direct solvers, where q is the bandwidth of the matrix and n is the degrees of freedom. The Fourier method is O ( N log( N )) in time when using the Fast Fourier Transform. We will perform an empirical evaluation of the e ffi ciency gained from the FM-SEM coupling compared to running the SEM in the entire domain for t max = 0 . 2 s. The methods are implemented

Figure 5: h-Convergence plots: FM{2,4,8,16,32}-SEM{4,8,16,32,64} are plotted pair-wise in the two top-most graphs, where 1) the temporal resolution is fixed for the blue graph, and 2) the temporal SEM resolution is twice the FM resolution for the blue graph. The lower red graph 3) shows convergence according to polynomial order P = 4 when oversampling 8 times the SEM P = 1 interface layer compared to the main SEM P = 4.

FM-SEM P=4 convergence

10 -1

10 -2

10 -3

L 2 error

10 -4

10 -5

10 -6

10 -2 10 -1

x [m]

Table 2: FM4-SEM8 CPU timings for l = 50 m

FM partition size (r) FM solver SEM solver interface

50% 1.2% 84.9% 13.0%

80% 6.0% 53.6% 39.6%

95% 10.1% 18.5% 70.3%

in Matlab, and the timings exclude the matrix assembly. The vast majority of the time in the SEM is spent solving the linear system of equations and is implemented using the Matlab backslash operator x = A \ b for solving the linear system Ax = b . For the FM, most of the work is spent in the Fourier transformation, where the build-in Matlab function ‘FFT’ is used. In all experiments, we will compare against the baseline forth-order SEM with Neumann boundaries. Future work should investigate more optimal solvers – such as sparse solvers – for the SEM to make the comparison more fair. In the first experiment, we compare the CPU time for a 1D case separately for the FM and SEM running in the entire 1D domain of lengths l = [6 , 13 , 25 , 50] m, i.e., with no couplings. The result is shown in Figure 6a, and we see 3x to 71x speedups depending on the domain size. The CPU times scale with O ( l 2 ) for the SEM and below the theoretical O ( l log( l )) limit for the FM. In the second experiment, we fix the 1D domain size to 50 m and compare the CPU time for the FM-SEM coupling for di ff erent Fourier partition sizes of r = [10 , 20 , 50 , 70 , 80 , 90 , 95] % relative to the entire domain. In Figure 6b the results are depicted, and we achieve speedups between 2x and 17x for Fourier partition sizes above 50% with scaling close to O ( r 3 ). In Table 2, the timings for l = 50 m are shown separately for the FM and SEM solvers and the interface handling. We notice that the SEM workload for a 50 / 50 partition split is taking 85% of the total computation time. Increasing the partition size of FM drastically decreases the SEM workload, and for r = 95 %, the workload of the interface handling starts dominating, taking up 70% of the time. Most of the workload at the interface is because of the space and time spline interpolation, though there is significant overhead when calling the interpolation methods interpolating only a few points near the interface. In fact, interpolating all pressure values instead of only the values around the interface has a minor impact on the absolute performance. Therefore, we expect more time-e ffi cient interpolations when extending to 2D and 3D, where much bigger pressure grids are to be interpolated near the interface.

Figure 6: CPU timings in a 1D domain. a) Comparison between the FM and the SEM running in the full domain for di ff erent domain sizes (no couplings), b) Comparison between FM-SEM and a baseline SEM with di ff erent relative FM partition sizes.

CPU time: FM vs. SEM in full domain

CPU time: FM-SEM coupling in L=50 m domain

10 3

0.8x

1.0x

10 2

10 2

2.8x

Time [s]

Time [s]

7.1x

10 1

9.7x

13.9x

3x 7x 21x 71x

10 1

17.8x

10 0

10 15 20 25 30 35 40 45 50 Length (l) [m]

10 20 30 40 50 60 70 80 90 FM partition size (r) [%]

(a)

(b)

CONCLUSION

We have implemented and coupled SEM and FM using a (2,6) FDTD interfacing scheme handling independent spatial and temporal resolutions. Coupling the Fourier method using four points per wavelength with the SEM using eight points per wavelength in a 5 m + 5m domain results in -36 dB interface errors and 9.1% relative errors compared to a reference solution. Using 12 points per wavelength for the SEM gives slightly better relative errors of 5.7%. The e ffi ciency of the coupled method is compared against an SEM running in the entire domain. For a fixed domain size of 50 m, the e ffi ciency of the coupled method is compared for di ff erent relative FM partition sizes, with an 18 times performance gain when 95% of the domain is running the FM. A more significant performance gain is expected to be achieved when going to larger 2D and 3D domains due to more degrees of freedom to be handled by the SEM. However, the workload at the interface will also grow compared to 1D, but we expect it to still be negligible compared to the saved computation time running the more expensive SEM solver. Future work could consider the use of more accurate time-stepping schemes such as the Runge- Kutta methods, improving the numerical accuracy of the SEM method. A current limitation of the implemented method is the need for an additional SEM layer of first-order polynomials dominating the overall convergence rate.

REFERENCES

[1] N. H. Asmar. Partial Di ff erential equations with Fourier Series and Boundary Value Problems . Pearson, Prentice Hall, 2 edition, 2005. [2] C. De Boor. A practical guide to splines , volume 27. Springer, 2001. [3] V. Dolean, P. Jolivet, and F. Nataf. An introduction to domain decomposition methods . Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2015. Algorithms, theory, and parallel implementation. [4] A. P. Engsig-Karup, C. Eskilsson, and D. Bigoni. A stabilised nodal spectral element method for fully nonlinear water waves. Journal of Computational Physics , 318:1–21, 2016.

[5] B. Hamilton and S. Bilbao. FDTD Methods for 3-D Room Acoustics Simulation with High- Order Accuracy in Space and Time. IEEE / ACM Transactions on Audio Speech and Language Processing , 25(11):2112–2124, 2017. [6] J. A. Hargreaves, L. R. Rendell, and Y. W. Lam. A framework for auralization of boundary element method simulations including source and receiver directivity. The Journal of the Acoustical Society of America , 145(4):2625–2637, 2019. [7] J. S. Hesthaven and T. Warburton. Nodal discontinuous Galerkin methods : algorithms, analysis, and applications , volume 54. Springer, 2008. [8] M. Hornikx, R. Waxler, and J. Forssén. The extended Fourier pseudospectral time-domain method for atmospheric sound propagation. The Journal of the Acoustical Society of America , 128(4):1632–1646, 2010. [9] H. Kuttru ff . Room Acoustics . CRC Press, 6 edition, 2016. [10] R. Mehra, N. Raghuvanshi, L. Savioja, M. C. Lin, and D. Manocha. An e ffi cient GPU-based time domain solver for the acoustic wave equation. Applied Acoustics , 73(2):83–94, 2012. [11] A. Melander, E. Strøm, F. Pind, A. P. Engsig-Karup, C.-H. Jeong, T. Warburton, N. Chalmers, and J. S. Hesthaven. Massive parallel nodal discontinuous galerkin finite element method simulator for room acoustics. Preprint on EPFL server , 2020. [12] T. Okuzono, T. Otsuru, R. Tomiku, and N. Okamoto. A finite-element method using dispersion reduced spline elements for room acoustics simulation. Applied Acoustics , 79:1–8, 2014. [13] R. Pagán Muñoz and M. Hornikx. Hybrid Fourier pseudospectral / discontinuous Galerkin time- domain method for wave propagation. Journal of Computational Physics , 348:416–432, 2017. [14] F. Pind, A. P. Engsig-Karup, C.-H. Jeong, J. S. Hesthaven, M. S. Mejling, and J. Strømann- Andersen. Time domain room acoustic simulations using the spectral element method. The Journal of the Acoustical Society of America , 145(6):3299–3310, 2019. [15] N. Raghuvanshi, N. Galoppo, and M. C. Lin. Accelerated wave-based acoustics simulation. In SPM ’08: Proceedings of the 2008 ACM symposium on Solid and physical modeling , pages 91–102. ACM, 2008. [16] N. Raghuvanshi, R. Narain, and M. C. Lin. E ffi cient and accurate sound propagation using adaptive rectangular decomposition. IEEE Transactions on Visualization and Computer Graphics , 15:789–801, September 2009. [17] B. Van Genechten, O. Atak, B. Bergen, E. Deckers, S. Jonckheere, J. S. Lee, A. Maressa, K. Vergote, B. Pluymers, D. Vandepitte, and W. Desmet. An e ffi cient Wave Based Method for solving Helmholtz problems in three-dimensional bounded domains. Engineering Analysis with Boundary Elements , 36(1):63–75, 2012. [18] B. Van Hal, W. Desmet, D. Vandepitte, and P. Sas. Hybrid Finite Element - Wave Based Method for acoustic problems. Computer Assisted Mechanics and Engineering Sciences , 10(4):479–494, 2003. [19] H. Xu, C. D. Cantwell, C. Monteserin, C. Eskilsson, A. P. Engsig-Karup, and S. J. Sherwin. Spectral / hp element methods: Recent developments, applications, and perspectives. Journal of Hydrodynamics , 30(1):1–22, 2018.