A Parallel Dissipation-Free and Dispersion-Optimized Explicit Time-Domain FEM for Large-Scale Room Acoustics Simulation

Abstract

Wave-based acoustics simulation methods such as finite element method (FEM) are reliable computer simulation tools for predicting acoustics in architectural spaces. Nevertheless, their application to practical room acoustics design is difficult because of their high computational costs. Therefore, we propose herein a parallel wave-based acoustics simulation method using dissipation-free and dispersion-optimized explicit time-domain FEM (TD-FEM) for simulating room acoustics at large-scale scenes. It can model sound absorbers with locally reacting frequency-dependent impedance boundary conditions (BCs). The method can use domain decomposition method (DDM)-based parallel computing to compute acoustics in large rooms at kilohertz frequencies. After validation studies of the proposed method via impedance tube and small cubic room problems including frequency-dependent impedance BCs of two porous type sound absorbers and a Helmholtz type sound absorber, the efficiency of the method against two implicit TD-FEMs was assessed. Faster computations and equivalent accuracy were achieved. Finally, acoustics simulation of an auditorium of 2271 ${\mathrm{m}}^3$ presenting a problem size of about 150,000,000 degrees of freedom demonstrated the practicality of the DDM-based parallel solver. Using 512 CPU cores on a parallel computer system, the proposed parallel solver can compute impulse responses with 3 s time length, including frequency components up to 3 kHz within 9000 s.

1. Introduction

1.1. Background

Wave-based acoustics simulation methods can model acoustics in architectural spaces [1] such as concert halls and classrooms accurately based on the wave equation or the Helmholtz equation using numerical techniques such as the finite element method (FEM) [2], finite difference method [3,4], and boundary element method [5]. A benefit of wave-based simulations is that they can accurately assess effects of two key technologies of sound absorption and sound scattering for room acoustics control [6]. With sound absorbers and acoustic diffusers, one can design a comfortable indoor sound environment by controlling noise levels and reverberation times according to room use. Consequently, wave-based simulations are reliable prediction tools for room acoustics design. Additionally, such simulations can further accommodate virtual reality (VR) technology. The use of such simulations in VR for engineering [7,8] and entertainment [9,10] is becoming an attractive application in recent years. An attractive feature of VR with wave-based simulations is that spatial sound field reproduction by ambisonics [11,12] with binaural rendering are conducted easily. That is true because expansion coefficients for spherical harmonics are accessed by placement of virtual microphone array [13] or numerical differentiation [14]. Subsequently, decoding and binaural auralization with headphones are accomplished using virtual studio technology plugins such as Resonance Audio [15] and SPARTA [16].

However, the main difficulty in applying wave-based simulations to architectural acoustic simulations is the necessary imposition of expensive computational loads, i.e., colossal memory consumption and long computational times, because architectural spaces have large dimensions compared to the acoustic wavelength at audible frequencies. Wave-based simulations require discretization of spaces with smaller-sized elements than the acoustic wavelength of interest to keep the discretization error, designated as the dispersion error, within an acceptable level. Although continuous progress in computing technology widens the applicable range of wave-based simulations, more efficient wave-based simulation methods able to reduce the dispersion error with smaller computational efforts than existing methods are necessary for application to architectural acoustic design. To this end, many researchers have examined the development of efficient wave-based simulation methods to date. For example, one can find reports of studies using higher-order finite-difference time-domain (FDTD) methods [17,18,19] and higher-order frequency-domain and time-domain FEMs [20,21,22,23,24,25,26]. Among other methods, this study addresses a time-domain FEM having explicit time-marching scheme, with the proposal of efficient acoustic simulation methods for large-scale room acoustics simulations.

Some efficient wave-based room acoustics simulation methods with dispersion-reduced FEMs using low-order finite elements (FEs) in both the frequency domain and time domain have been developed [27,28,29,30,31,32,33,34]. They have been applied respectively to acoustics simulation in a simple-shaped concert hall [31], simulations of reverberation absorption coefficient measurements of sound absorbers with woven and non-woven fabrics [32], and acoustics improvement of an existing meeting room with microperforated panel absorbers [33]. The salient feature of the dispersion-reduced FEMs is that they can maintain dispersion error at a low level with low-order FEs without increasing additional computational costs. The methods are able to do so because of the use of modified integration rule (MIR) [35,36], which uses numerical integration points able to reduce dispersion error for element matrix construction. Consequently, the methods can decrease the total number of degrees of freedom (DOFs) without increasing the matrix’s bandwidth more than higher-order FEM. The computational costs required per element are lower than those of higher-order methods. For development, first, dispersion-reduced implicit time-domain FEM (TD-FEM) [27] and frequency-domain FEM (FD-FEM) [29] have been developed for room acoustics simulation. A dispersion-reduced high-order FD-FEM and TD-FEM using spline FEs have also been proposed in an earlier report of the relevant literature [26]. Subsequently, a dispersion-reduced explicit TD-FEM with fourth-order accuracy in terms of dispersion error was proposed to accelerate a computational speed for acoustic simulations further in architectural spaces discretized only with rectangular FEs [28,30]. Later, by developing a time integration method called the modified Adams method for time discretization [34], the limitation of using rectangular FEs was removed without accuracy reduction caused by dispersion error. However, the explicit TD-FEM with the modified Adams method still has the important shortcoming that includes numerical dissipation error. With the error, sound energy decreases over time despite the use of a lossless wave equation. Although the attenuation level in sound energy can be reduced using a smaller time interval value, the strategy entails increased computational costs.

More recently, a dissipation-free fourth-order accurate explicit TD-FEM for room acoustics simulation [37] has been proposed, which overcomes shortcomings in TD-FEM with the modified Adams method. It achieves a dissipation-free scheme using a newly developed three-step time-integration method based on the general form of linear multi-step time integration method [38]. Furthermore, the dissipation-free explicit TD-FEM can enhance the capability of sound field approximation at higher frequencies using an optimization method with no additional computational costs. In the optimized method, the dispersion errors in axial and diagonal directions are minimized at a specific mesh or element resolution with optimized MIR and weight coefficients in time integration. The basic performance of the dissipation-free and dispersion-optimized explicit TD-FEMs has been examined for simply shaped room acoustics problems with rigid boundaries at wideband frequencies. Those examinations showed outstanding performance against a standard second-order accurate implicit TD-FEM. Nevertheless, the method still lacks implementation of absorbing boundary conditions (BCs) to reflect the effects of sound absorbers such as porous materials. Therefore, its performance for practical room acoustic modeling remains unclear.

The quality of sound absorption modeling is well known to have a strong effect on the resulting accuracy in sound field predictions. Ideally, using an extended-reaction BC ability is desirable to address both frequency dependence and incident-angle dependence in sound absorption characteristics of sound absorbers. Note that uncertainty in input data concerning sound absorption modeling also has a considerable effect on calculation results [39,40]. A local-reaction BC that simplifies the incident-angle dependence is another selection because the computation becomes easier with the surface impedance of sound absorbers than the extended-reaction BCs that require discretization of materials according to corresponding partial differential equations. Although the implementation of frequency-dependence into time-domain methods is difficult because of the inclusion of convolution operation, very recently developed time-domain methods [23,41,42,43,44,45] can specifically address the frequency-dependence in the local-reaction BC and extended-reaction BC efficiently. This report also addresses the implementation of frequency-dependent local-reaction impedance BCs into the dissipation-free and dispersion-optimized explicit TD-FEM using auxiliary differential equations method [46].

Furthermore, the use of parallel computing technology for wave-based room acoustics simulation is a simple choice for extending its applicable range to large architectural spaces at high frequencies. Because cloud-computing environments such as Amazon Web Service [47], Microsoft Azure [48], and Google Cloud [49] are becoming increasingly familiar in recent years, the development of massively parallel solvers is an exciting subject. Some earlier studies have presented excellent parallel performance of parallel wave-based acoustics simulation with a high performance computing (HPC) system having many CPU cores or GPUs [50,51,52,53,54,55]. A marked speeding up of acoustics simulation with a parallel implicit TD-FEM using spline FEs have also been demonstrated for a simply shaped hall with 13,000 ${\mathrm{m}}^3$ [50] when using an HPC system with 256 cores about ten years ago. Because of this marked performance gain with large-scale parallel computing can also anticipate for a dissipation-free and dispersion-optimized explicit TD-FEM, as a practical demonstration, we will examine the parallel performance for calculating an impulse response in an auditorium model, including frequency-dependent impedance BCs of two porous-type sound absorbers and a Helmholtz-type sound absorber.

1.2. Purpose of This Study

The present study proposes a parallel wave-based room-acoustics simulation method with a dissipation-free and dispersion-optimized explicit TD-FEM, with eventual demonstration of its practicality for acoustics simulation at kilohertz frequency range in a large room. As the salient point of novelty of the present work, we present the completely new simulation method for large-scale room acoustics simulation composed of a time-marching scheme to deal with frequency-dependent impedance BC and a DDM-based parallel computation algorithm. It can predict an impulse response up to three seconds in an auditorium of 2271 ${\mathrm{m}}^3$ having multiple frequency-dependent impedance BCs within 9000 s using 512 cores. Our study also includes two validation studies that use impedance tube problems and a small cubic room problem with frequency-dependent impedance BCs. For the latter problem, we show the performance of the proposed method over existing TD-FEMs and the efficient numerical setup for the proposed method. Hereinafter, the proposed dissipation-free and dispersion-optimized explicit TD-FEM and the non-optimized fourth-order explicit TD-FEM, i.e., the dissipation-free method only, are abbreviated respectively to the Opt-E TD-FEM and the 4th-E TD-FEM. This paper is organized as follows. Section 2 describes the theory of proposed time-marching scheme with Opt-E and 4th-E TD-FEMs including frequency-dependent impedance BCs. This section also includes, for readers’ convenience, details of theories of other FEM-based room acoustic solvers used herein. Section 3 presents tests of the validity of the proposed methods with the impedance tube problem and cubic room problems. Section 4 examines the efficiency of the proposed method over two existing implicit TD-FEM, including a proposal of the efficient optimization setup for the proposed methods. Section 5 demonstrates the practicality of parallel Opt-E TD-FEM using a hybrid framework with Message Passing Interface (MPI) and OpenMP to large-scale room acoustics simulation with an auditorium model with the volume of 2271 ${\mathrm{m}}^3$. Section 6 concludes the paper. All computational codes used in the present study are in-house software written with Fortran 90.

2. Theory

This section presents the theory of proposed room acoustics simulation methods with the dissipation-free and dispersion-optimized explicit TD-FEMs: Opt-E and 4th-E TD-FEMs. For the readers’ convenience, this section also includes brief descriptions of the fundamental theories of FD-FEM and implicit TD-FEM used for subsequent numerical experiments. The end of this section includes presentation of the theoretical dispersion error properties of the proposed Opt-E and 4th-E TD-FEMs compared to two implicit TD-FEMs. The results are expected to be useful for discussion in subsequent numerical experiments.

2.1. Basic Equation and Its Spatial Discretization with Finite Elements

To model acoustic wave propagation in a space $\Omega$ enclosed by a boundary $\Gamma$, we specifically examine the second-order scalar wave equation expressed as

$$\frac{{\partial }^{2}p(\mathit{r},t)}{\partial t^2}-c_0^2{\nabla }^{2}p(\mathit{r},t)={\rho }_0{c}_0^2\frac{\partial q(\mathit{r},t)}{\partial t}\delta (\mathit{r}-{\mathit{r}}_{\mathrm{s}})$$

(1)

with the following three BCs: the rigid BC ${\Gamma }_0$, the vibration BC ${\Gamma }_{\mathrm{V}}$, and the impedance BC ${\Gamma }_{\mathrm{Z}}$ as

$$\frac{\partial p(\mathit{r},t)}{\partial n}=\left\{\begin{array}{cc}0\hfill & \mathrm{on}{\Gamma }_0\hfill \\ -{\rho }_0{\dot{v}}_{\mathrm{n}}(\mathit{r},t)\hfill & \mathrm{on}{\Gamma }_{\mathrm{V}}\hfill \\ -\frac{1}{c_0}{\int }_{-\infty }^t{y}_{\mathrm{n}}(\mathit{r},t-\tau )\dot{p}(\mathit{r},\tau )\mathrm{d}\tau \hfill & \mathrm{on}{\Gamma }_{\mathrm{Z}}\hfill \end{array}\right.$$

(2)

Therein, p stands for the sound pressure, q signifies the volume velocity per unit volume, and t denotes the time. Also, ${\nabla }^2$ is the Laplacian operator. $c_0$ represents the speed of sound. ${\rho }_0$ denotes the density of air. For the analyses described herein, we consistently set $c_0$ = 343.7 m/s and ${\rho }_0$ = 1.205 kg/m${}^3$. Also, $\mathit{r}$ and ${\mathit{r}}_s$ respectively represent the arbitrary coordinate vector in $\Omega$ and the coordinate vector of a point source. The Dirac delta function is $\delta$. In addition, $v_{\mathrm{n}}$ represents the vibration velocity. Symbol · is the first-order time derivative. $y_n$ stands for the specific acoustic admittance ratio in the time domain. Multiplication of Equation (1) by the arbitrary weight function ${\varphi }_{\mathrm{f}}$ and applying Green’s theorem yields the weak form defined as

$${\int }_{\Omega }{\varphi }_{\mathrm{f}}\frac{{\partial }^{2}p(\mathit{r},t)}{\partial t^2}+c_0^2\nabla {\varphi }_{\mathrm{f}}\nabla p(\mathit{r},t)\mathrm{d}\Omega -c_0^2{\int }_{\Gamma }{\varphi }_{\mathrm{f}}\frac{\partial p(\mathit{r},t)}{\partial n}\mathrm{d}\Gamma ={\rho }_0{c}_0^2{\int }_{\Omega }{\varphi }_{\mathrm{f}}\frac{\partial q(\mathit{r},t)}{\partial t}\delta (\mathit{r}-{\mathit{r}}_{\mathrm{s}})\mathrm{d}\Omega$$

(3)

Taking the Galerkin method with three BCs in Equation (2) leads to the second-order ordinary differential equation (ODE) as [45]

$$\mathit{M}\ddot{\mathit{p}}+c_0^2\mathit{K}\mathit{p}+c_0\mathit{C}(\check{y}\times \dot{\mathit{p}})=\mathit{f}$$

(4)

where

$$\begin{array}{ccc}& & \mathit{M}=\sum _{i=1}^{N_{\mathrm{e}}}{\mathit{M}}_{\mathrm{e}}=\sum _{i=1}^{N_{\mathrm{e}}}{\int }_{{\Omega }_{\mathrm{e}}}{\mathit{N}}^T\mathit{N}\mathrm{d}{\Omega }_{\mathrm{e}}\hfill \end{array}$$

(5)

$$\begin{array}{ccc}& & \mathit{K}=\sum _{i=1}^{N_{\mathrm{e}}}{\mathit{K}}_{\mathrm{e}}=\sum _{i=1}^{N_{\mathrm{e}}}{\int }_{{\Omega }_{\mathrm{e}}}\nabla {\mathit{N}}^T\nabla \mathit{N}\mathrm{d}{\Omega }_{\mathrm{e}}\hfill \end{array}$$

(6)

$$\begin{array}{ccc}& & \mathit{C}=\sum _{i=1}^{N_{\mathrm{s}}}{\mathit{C}}_{\mathrm{e}}=\sum _{i=1}^{N_{\mathrm{s}}}{\int }_{{\Gamma }_{\mathrm{e}}}{\mathit{N}}^T\mathit{N}\mathrm{d}{\Gamma }_{\mathrm{e}}\hfill \end{array}$$

(7)

$$\begin{array}{ccc}& & \check{y}\times \dot{\mathit{p}}={\int }_{-\infty }^t{y}_n(\mathit{r},t-\tau )\dot{\mathit{p}}(\mathit{r},\tau )\mathrm{d}\tau \hfill \end{array}$$

(8)

Here, $\mathit{M}$, $\mathit{K}$, and $\mathit{C}$ respectively represent the global consistent mass matrix, the global stiffness matrix, and the global dissipation matrix composed with each element matrix ${\mathit{M}}_{\mathrm{e}}$, ${\mathit{K}}_{\mathrm{e}}$ and ${\mathit{C}}_{\mathrm{e}}$. Also, $\mathit{N}$ denotes the shape function vector; $\mathit{p}$ and $\mathit{f}$ respectively denote the sound pressure vector and the external force vector. $N_{\mathrm{e}}$ represents the number of FE ${\Omega }_{\mathrm{e}}$; $N_{\mathrm{s}}$ is the number of surface elements ${\Gamma }_{\mathrm{e}}$. Symbol $··$ signifies second-order time derivatives. For the analyses described herein, only eight-node hexahedral FEs are used for spatial discretization. Therefore, integrations for element matrix calculations in Equations (5) and (6) are performed with two-point Gauss–Legendre quadrature in each direction as

$$\begin{array}{ccc}& & {\mathit{K}}_{\mathrm{e}}=\sum _{h=1}^2\sum _{i=1}^2\sum _{j=1}^2\nabla \mathit{N}{({\alpha }_{k,h},{\alpha }_{k,i},{\alpha }_{k,j})}^{\mathrm{T}}\nabla \mathit{N}({\alpha }_{k,h},{\alpha }_{k,i},{\alpha }_{k,j})det(\mathit{J})\hfill \end{array}$$

(9)

$$\begin{array}{ccc}& & {\mathit{M}}_{\mathrm{e}}=\sum _{h=1}^2\sum _{i=1}^2\sum _{j=1}^2\mathit{N}{({\alpha }_{m,h},{\alpha }_{m,i},{\alpha }_{m,j})}^{\mathrm{T}}\mathit{N}({\alpha }_{m,h},{\alpha }_{m,i},{\alpha }_{m,j})det(\mathit{J})\hfill \end{array}$$

(10)

where ${\alpha }_k$ and ${\alpha }_m$ represent coordinates of integration points in a natural coordinate system for ${\mathit{K}}_{\mathrm{e}}$ and ${\mathit{M}}_{\mathrm{e}}$. Also, $\mathit{J}$ is the Jacobian matrix. The standard acoustic FEMs in both frequency-domain and time-domain present second-order accuracy in space with conventional integration point coordinates of ${\alpha }_k={\alpha }_m=\pm \sqrt{1/3}$. However, it can achieve fourth-order accuracy in space using modified coordinates of ${\alpha }_k={\alpha }_m=\pm \sqrt{2/3}$ [27,31,35,36]. Applying the Fourier transformation to Equation (4) with the time factor of $exp(\mathrm{j}\omega t)$ yields the following discretized matrix equation for FD-FEM.

$$\left(\mathit{K}-k^2\mathit{M}+\mathrm{j}k\mathit{C}{y}_{\mathrm{n}}(\omega )\right)\mathit{p}=\mathit{f}$$

(11)

Therein, k and $\mathrm{j}$ respectively denote the wavenumber and the imaginary unit. Furthermore, $\omega$ represents the angular frequency and $y_{\mathrm{n}}(\omega )$ denotes the frequency-domain specific acoustic admittance ratio. The FD-FEM can simulate sound fields at a steady state by solving Equation (11) with sparse direct solvers such as PARDISO or Krylov-subspace iterative solvers such as Conjugate Orthogonal Conjugate Gradient (COCG) method [56]. In later numerical experiments Equation (11) are solved using the COCG method with a diagonal scaling preconditioning under the convergence tolerance of 10${}^{-6}$.

2.2. Incorporation of Frequency-Dependent Absorbing Boundary

For the analyses presented herein, the ADE method [46] is used for dealing with convolution involved with the frequency-dependent impedance boundary ${\Gamma }_{\mathrm{Z}}$ in Equation (2) efficiently. The ADE method uses the frequency-domain specific acoustic admittance ratio $y_{\mathrm{n}}(\omega )$, which is approximated to a rational form as

$$y_{\mathrm{n}}(\omega )\approx y^{\prime }(\omega )=y_{\infty }+\sum _{i=1}^{n_{\mathrm{rp}}}\frac{A_i}{{\lambda }_i+\mathrm{j}\omega }+\sum _{i=1}^{n_{\mathrm{cp}}}(\frac{B_i-\mathrm{j}{C}_i}{{\alpha }_i-\mathrm{j}{\beta }_i+j\omega }+\frac{B_i+\mathrm{j}{C}_i}{{\alpha }_i+\mathrm{j}{\beta }_i+\mathrm{j}\omega })$$

(12)

where $n_{\mathrm{rp}}$ and $n_{\mathrm{cp}}$ are the numbers of real poles ${\lambda }_i$ and pairs of complex conjugate poles ${\alpha }_i\pm \mathrm{j}{\beta }_i$. Also, $y_{\infty }$, $A_i$, $B_i$ and $C_i$ are real-valued parameters. To perform stable simulations ${\lambda }_i\ge 0$, ${\alpha }_i\ge 0$ are imposed for causality reasons. Also, $\mathrm{Re}\left(y^{\prime }(\omega )\right)\ge 0$ at arbitrary $\omega$ is imposed for passivity. The use of $y^{\prime }(\omega )$ approximates the convolution in Equation (4) as

$$\check{y}\ast \dot{\mathit{p}}\approx y_{\infty }\dot{\mathit{p}}+\sum _{i=1}^{n_{\mathrm{rp}}}{A}_i{\mathit{\varphi }}_i+2\sum _{i=1}^{n_{\mathrm{cp}}}(B_i{\mathit{\psi }}_i^{(1)}+C_i{\mathit{\psi }}_i^{(2)})$$

(13)

with accumulators ${\mathit{\varphi }}_i$, ${\mathit{\psi }}_i^{(1)}$ and ${\mathit{\psi }}_i^{(2)}$ respectively denoted by the following simultaneous first-order ODEs as

$$\begin{array}{ccc}& & {\dot{\mathit{\varphi }}}_i+{\lambda }_i{\mathit{\varphi }}_i=\dot{\mathit{p}}\hfill \end{array}$$

(14)

$$\begin{array}{ccc}& & {\dot{\mathit{\psi }}}_i^{(1)}+{\alpha }_i{\mathit{\psi }}_i^{(1)}+{\beta }_i{\mathit{\psi }}_i^{(2)}=\dot{\mathit{p}}\hfill \end{array}$$

(15)

$$\begin{array}{ccc}& & {\dot{\mathit{\psi }}}_i^{(2)}+{\alpha }_i{\mathit{\psi }}_i^{(2)}-{\beta }_i{\mathit{\psi }}_i^{(1)}=\mathbf{0}\hfill \end{array}$$

(16)

This report presents a solution for the ODEs of Equations (14), (15) and (16) above using the Crank–Nicolson method in both implicit and explicit TD-FEMs. As shown in numerical simulations conducted in Section 3, the Crank–Nicolson method has good accuracy and stability in the time-marching calculations for the accumulators. Details of the Crank–Nicolson method implementation are presented in an earlier report of the literature [45]. With the help of ADE method, Equation (4) is transformed into

$$\mathit{M}\ddot{\mathit{p}}+c_0^2\mathit{K}\mathit{p}+c_0\mathit{C}(y_{\infty }\dot{\mathit{p}}+\sum _{i=1}^{n_{\mathrm{rp}}}{A}_i{\mathit{\varphi }}_i+2\sum _{i=1}^{n_{\mathrm{cp}}}(B_i{\mathit{\psi }}_i^{(1)}+C_i{\mathit{\psi }}_i^{(2)}))=\mathit{f}.$$

(17)

The implicit TD-FEM solves Equation (17) with a direct time integration method. In a later numerical experiment, fourth-order accurate Fox–Goodwin method [57] is used for time integration. The stability limit time interval, $\Delta t_{\mathrm{limit}}$ of the implicit TD-FEM with Fox–Goodwin method depends on the local coordinates of integration points as shown below [27,31].

$$\begin{array}{ccc}& & \Delta t_{\mathrm{limit}}=\frac{h}{\sqrt{6}{c}_0}\mathrm{with}{\alpha }_k={\alpha }_m=\pm \sqrt{\frac{1}{3}}\hfill \end{array}$$

(18)

$$\begin{array}{ccc}& & \Delta t_{\mathrm{limit}}=\frac{h}{\sqrt{3}{c}_0}\mathrm{with}{\alpha }_k={\alpha }_m=\pm \sqrt{\frac{2}{3}}\hfill \end{array}$$

(19)

In those equations, h denotes the edge size of cubic-shaped FE. Hereinafter, we call the implicit TD-FEM with integration points in Equation (18) 2nd-I TD-FEM, and those with integration points in Equation (19) 4th-I TD-FEM. The 2nd-I TD-FEM has second-order accuracy in space and fourth-order accuracy in time, whereas the 4th-I TD-FEM achieves fourth-order accuracy in both space and time. The two methods have the same dispersion error property with respect to temporal discretization, but the 2nd-I TD-FEM has larger dispersion error than the 4th-I TD-FEM because dispersion error derives from spatial discretization, as shown in Section 2.4. To solve the linear system of equations at each time step, we use an iterative solver, conjugate gradient (CG) method, with a diagonal scaling preconditioning for both implicit TD-FEMs. Its convergence tolerance is set as 10${}^{-6}$.

2.3. Proposed Dissipation-Free and Dispersion-Optimized Explicit TD-FEMs including Frequency-Dependent Impedance Boundary

2.3.1. Dissipation-Free Time-Stepping Scheme

The explicit TD-FEM solves the following simultaneous first-order ODEs equivalent to the second-order system of Equation (4) [34].

$$\begin{array}{ccc}& & \mathit{D}\dot{\mathit{p}}=\mathit{M}\mathit{v}\hfill \end{array}$$

(20)

$$\begin{array}{ccc}& & \mathit{D}\dot{\mathit{v}}=\mathit{f}-c_0^2\mathit{K}\mathit{p}-c_0\mathit{C}(y_{\infty }\dot{\mathit{p}}+\sum _{i=1}^{n_{\mathrm{rp}}}{A}_i{\mathit{\varphi }}_i+2\sum _{i=1}^{n_{\mathrm{cp}}}(B_i{\mathit{\psi }}_i^{(1)}+C_i{\mathit{\psi }}_i^{(2)}))\hfill \end{array}$$

(21)

In those equations, $\mathit{D}$ represents the lumped mass matrix assembled from $\mathit{M}$ using the row-sum method [58]. Additionally, $\mathit{v}$ denotes the auxiliary vector equivalent to the first derivative of $\mathit{p}$ with respect to time. To achieve higher accuracy in time without numerical dissipation, we discretize $\dot{\mathit{p}}$ in Equation (20) using the three-step time integration method [37], and $\dot{\mathit{v}}$ in Equation (21) using first-order backward difference approximation as

$$\begin{array}{ccc}& & {\mathit{p}}^n=\sum _{i=1}^3{a}_i{\mathit{p}}^{n-i}+\Delta t{\mathit{D}}^{-1}\sum _{j=1}^3{b}_j{\mathit{Mv}}^{n-j}\hfill \end{array}$$

(22)

$$\begin{array}{ccc}\hfill & & {\mathit{v}}^n={\mathit{v}}^{n-1}+\Delta t{\mathit{D}}^{-1}\left({\mathit{f}}^n-c_0^2\mathit{K}{\mathit{p}}^n-c_0\mathit{C}(y_{\infty }\dot{\mathit{p}}+\sum _{i=1}^{n_{\mathrm{rp}}}{A}_i{\mathit{\varphi }}_i+2\sum _{i=1}^{n_{\mathrm{cp}}}(B_i{\mathit{\psi }}_i^{(1)}+C_i{\mathit{\psi }}_i^{(2)}))\right)\hfill \end{array}$$

(23)

Here, $\Delta t$ and n respectively denote the time interval and the time step. Furthermore, $a_i$ and $b_j$ respectively represent the i-th and j-th real-valued weight coefficients. An inherently dissipation-free scheme is constructed with

$$a_1=2,a_2=-2,a_3=1,b_2=1-2b_1,b_3=b_1$$

(24)

Details of the mechanism involved in the elimination of dissipation error are described in an earlier report of the relevant literature [37]. Additionally, $\dot{\mathit{p}}$ in Equation (23) is discretized using first-order central difference approximation based on Equation (22) as

$$\begin{array}{cc}\hfill {\dot{\mathit{p}}}^n& =\frac{b_1}{2}{\mathit{D}}^{-1}{\mathit{Mv}}^n+\frac{1}{2}{\mathit{D}}^{-1}(b_2{\mathit{Mv}}^{n-1}+b_3{\mathit{Mv}}^{n-2})\hfill \\ \hfill & +\frac{1}{2\Delta t}(a_1{\mathit{p}}^n+(a_2-1){\mathit{p}}^{n-1}+a_3{\mathit{p}}^{n-2})\hfill \end{array}$$

(25)

Substituting Equation (25) into Equation (23) produces the time-marching equation for $\mathit{v}$ as

$$\begin{array}{cc}\hfill & \left(\mathit{I}+\frac{b_1}{2}{c}_0\Delta t{y}_{\infty }{\mathit{D}}^{-1}\mathit{C}{\mathit{D}}^{-1}\mathit{M}\right){\mathit{v}}^n\hfill \\ \hfill & ={\mathit{v}}^{n-1}+\Delta t{\mathit{D}}^{-1}\left({\mathit{f}}^n-c_0^2\mathit{K}{\mathit{p}}^n-c_0\mathit{C}(y_{\infty }{\mathit{V}}^n+\sum _{i=1}^{n_{\mathrm{rp}}}{A}_i{\mathit{\varphi }}_i+2\sum _{i=1}^{n_{\mathrm{cp}}}(B_i{\mathit{\psi }}_i^{(1)}+C_i{\mathit{\psi }}_i^{(2)}))\right)\hfill \end{array}$$

(26)

with

$${\mathit{V}}^n=\frac{1}{2}{\mathit{D}}^{-1}(b_2{\mathit{Mv}}^{n-1}+b_3{\mathit{Mv}}^{n-2})+\frac{1}{2\Delta t}(a_1{\mathit{p}}^n+(a_2-1){\mathit{p}}^{n-1}+a_3{\mathit{p}}^{n-2})$$

(27)

The dissipation-free time-stepping scheme for explicit TD-FEM is constructed with Equations (22) and (26). Consequently, both the 4th-E TD-FEM and the Opt-E TD-FEM use Equations (22) and (26) with the weight coefficients of Equation (24), but with different integration points of ${\alpha }_m$ and ${\alpha }_k$ and different weight coefficients of $b_1$. The 4th-E TD-FEM achieves fourth-order accuracy in space and time using ${\alpha }_m=\pm \sqrt{4/3}$, ${\alpha }_k=\pm \sqrt{2/3}$, and $b_1=13/12$. However, the Opt-E TD-FEM uses theoretically derived optimum integration points and weight coefficient, i.e., ${\alpha }_m=\pm {\alpha }_{m,\mathrm{opt}}$, ${\alpha }_k=\pm {\alpha }_{k,\mathrm{opt}}$ and $b_1=b_{1,\mathrm{opt}}$, respectively denoted as [37]

$${\alpha }_{m,\mathrm{opt}}=\sqrt{\frac{4{\pi }^2+R^2({cos}^2(\frac{2\pi }{R})-1)}{R^2{(cos(\frac{2\pi }{R})-1)}^2}}$$

(28)

$${\alpha }_{k,\mathrm{opt}}=\sqrt{\frac{-(\frac{2}{R}\sqrt{X_1}+X_2)}{3X_3^2{X}_4^2{X}_5{X}_6}}$$

(29)

with

$$X_1=-X_3^2{X}_4^2{X}_6(R^2{X}_4^2{(Y_1-Y_2)}^2(-X_6)+48{\pi }^2{X}_5)$$

$$\begin{array}{cc}\hfill X_2& ={\alpha }_m^6{X}_3^3{X}_5{X}_7-X_3{(Y_1+1)}^2(1+Y_2)X_7\hfill \\ \hfill & +{\alpha }_m^2(Y_1^2-1)(-3+Y_1-Y_2+3Y_1{Y}_2)X_7\hfill \\ \hfill & -{\alpha }_m^4{X}_3^2{X}_7^2\hfill \end{array}$$

$$\begin{array}{ccc}& & X_3=Y_1-1,X_4=1+{\alpha }_m^2+Y_1-{\alpha }_m^2{Y}_1\hfill \\ & & X_5=Y_2-1,X_6=1+{\alpha }_m^2+Y_2-{\alpha }_m^2{Y}_2\hfill \\ & & X_7=-3-Y_1-Y_2+3Y_1{Y}_2\hfill \end{array}$$

$$Y_1=cos\frac{\pi }{R},Y_2=cos\frac{\sqrt{2}\pi }{R}$$

(30)

and

$$b_{1,\mathrm{opt}}=\frac{1-2cos(\omega \Delta t)}{{\omega }^2\Delta t^2}+\frac{1}{4}{csc}^2(\frac{\omega \Delta t}{2})$$

(31)

In the equations presented above, R represents the spatial resolution of used mesh or each FE, which is expressed as the number of elements per wavelength at an arbitrary frequency. Therein ${\alpha }_{m,\mathrm{opt}}$ and ${\alpha }_{k,\mathrm{opt}}$ minimize the dispersion error at arbitrary frequencies under given spatial resolution meshes, respectively, in the axial direction and diagonal direction. In addition, $b_{1,\mathrm{opt}}$ minimizes the dispersion error in temporal discretization at the frequency of interest. Both 4th-E and Opt-E TD-FEMs have the same computational cost in terms of the required memory and computational complexity. However, as shown in Section 4, the Opt-E TD-FEM can enhance accuracy at higher frequencies compared to 4th-E TD-FEM. The stability conditions for both methods become the same; $\Delta t_{\mathrm{limit}}$ is given as

$$\Delta t_{\mathrm{limit}}=0.490774h/c_0$$

(32)

2.3.2. Locally Implicit Strategy

Here, Equation (22) is an explicit equation, but Equation (26) becomes an implicit equation, which requires the solution of linear systems because the consistent mass matrix $\mathit{M}$ of the left-hand side engenders non-diagonal components in the coefficient matrix. The fact vitally reduces the computational efficiency of explicit TD-FEM. However, the non-diagonal components in the coefficient matrix of Equation (26) appear only for row-related nodes on the boundary surfaces. Consequently, the local application of linear equation solver into the boundary part reduces the computational complexity accompanying the solution of linear equations. Using this strategy, the remaining part of the system is calculable explicitly using the following time-marching equation.

$${\mathit{v}}^n={\mathit{v}}^{n-1}+\Delta t{\mathit{D}}^{-1}\left({\mathit{f}}^n-c_0^2\mathit{K}{\mathit{p}}^n\right)$$

(33)

A report of earlier work [59] described that the locally implicit strategy can make the computational complexity much smaller than that in a case with a fully implicit scheme, but it can do so without reducing accuracy. We apply the locally implicit scheme for both the Opt-E and the 4th-E TD-FEMs. As a solver for local linear equations having an indefinite coefficient matrix, we used the conjugate residual (CR) method [60] without preconditioning. Subsequent numerical experiments reveal that CR shows rapid convergence with a small iteration number of less than ten.

2.4. Theoretical Dispersion Error Property of Proposed Explicit TD-FEMs against Two Implicit TD-FEMs

This section presents the theoretically estimated dispersion error property of the 4th-E and Opt-E TD-FEMs compared to those of the 4th-I and 2nd-I TD-FEMs. It is noteworthy that dispersion error estimation is based on an ideal condition that assumes plane wave propagation in free field on a structured mesh discretized using cubic FEs. The theoretical dispersion errors of the 4th-E method and two implicit TD-FEMs are evaluated explicitly with the equations written in [37] and [27], respectively. The dispersion error of Opt-E method can be evaluated numerically according to procedures described in our earlier work [37]. We evaluated the dispersion error of Opt-E method with four optimizing conditions, where each condition minimized the dispersion error at the respective spatial resolutions of R = 6.87, 6.25, 5.5, and 4.91. For all methods, time intervals $\Delta t$ were set to their stability limit values. Figure 1 presents dispersion errors of 4th-E and Opt-E TD-FEMs with the four optimized conditions compared to those of 4th-I and 2nd-I TD-FEMs. In the figure, positive error values numerically represent sound speeds lower than $c_0$, whereas negative values represent higher sound speeds. Additionally, the illustrated curves show errors in sounds propagating for Cartesian coordinate axes because all methods show maximum dispersion error in those directions.

The standard 2nd-I TD-FEM needs $R>12.8$ to maintain dispersion error within 1%. Other methods can suppress dispersion error of less than 1% with much fewer R: the 4th-E method requires $R>5.9$; the Opt-E method (optimized at $R=6.87$) requires $R>4.5$; the Opt-E method (optimized at $R=6.25$) needs $R>4.4$; the Opt-E method (optimized at $R=5.5$) needs $R>4.2$; the Opt-E method (optimized at $R=4.91$) needs $R>4.0$; and the 4th-I method requires $R>4.1$. We also observed that the Opt-E methods took minimum values at intended frequencies for the respective optimizing conditions. Moreover, the Opt-E TD-FEM keeps the dispersion error $e_{\mathrm{dispersion}}$ small in a wideband and at a higher frequency than 4th-E method, i.e., $e_{\mathrm{dispersion}}\le 0.2\%$ in $R\ge 5.8$ (optimized at R = 6.87), $e_{\mathrm{dispersion}}\le 0.3\%$ in $R\ge 5.3$ (optimized at R = 6.25), $e_{\mathrm{dispersion}}\le 0.4\%$ in $R\ge 5$ (optimized at R = 5.5), and $e_{\mathrm{dispersion}}\le 0.5\%$ in $R\ge 4.4$ (optimized at R = 4.91). In addition, the 4th-E method shows slower sound speed, but Opt-E methods show a faster sound speed at a spatial resolution larger than that used for optimization.

3. Validation of the Proposed Explicit TD-FEM with the Impedance Tube Problem and Small Cubic Room Problem

This section presents the validity of the proposed explicit TD-FEM, including frequency-dependent impedance BCs via numerical analyses of an impedance tube and a cubic room. The verification was made by comparison with a theoretical solution in the impedance tube problem and a reference solution accessed by the well-developed FD-FEM with a mesh having sufficient spatial resolution in the cubic room problem. Here, we use only the 4th-E TD-FEM, which is a basis of the Opt-E TD-FEM. The accuracy of Opt-E TD-FEM is examined in the next section. For the validation study, we modeled sound absorption characteristics of three sound absorbers for frequency-dependent impedance boundaries: rigid-backed glass-wool (GW), rigid-backed needle felt (NF), and glass-wool backed microperforated panel (MPP-GW). The acoustical characteristics of two porous absorbers were modeled using an equivalent fluid model with Miki model: GW and NF [61]. Also, GW was assumed to have flow resistivity of 6900 Pa s/m${}^2$ and 50 mm thickness, whereas NF has flow resistivity of 10,000 Pa s/m${}^2$ and 15 mm thickness. The surface impedance Z of each porous absorber was calculated using the transfer matrix method [62]. Additionally, we used Maa’s transfer impedance model [63] for MPP modeling. Then, Z for MPP-GW is also calculated based on the transfer matrix method. For the following experiments, we set the MPP’s parameters as the panel thickness of 0.5 mm, the hole diameter of 0.5 mm, the perforation ratio of 0.55% and the dynamic viscosity of the air of 17.9 $\mathsf{\mu }$Pa s.

The surface impedance of each absorber was used to calculate the parameters in the rational function form of Equation (12) using the vector fitting method [64] with passivity enforcement [65] at frequencies of 100 Hz to 5 kHz. The fittings were made respectively with $n_{\mathrm{rp}}$ = 3 and $n_{\mathrm{cp}}$ = 3 for GW, $n_{\mathrm{rp}}$ = 2 and $n_{\mathrm{cp}}$ = 2 for NF and $n_{\mathrm{rp}}$ = 10 and $n_{\mathrm{cp}}$ = 5 for MPP-GW. The complete parameters for rational function forms of each absorber are presented in Appendix A.

3.1. Impedance Tube Problem

This subsection shows the validity of 4th-E TD-FEM with an impedance tube problem. With 4th-E TD-FEM, we computed the normal incidence sound absorption characteristics of three sound absorbers, i.e., absorption coefficient and surface impedance, according to the transfer function method [66]. Then the computed results were compared with the theoretical results calculated using the transfer matrix method. Figure 2 presents the analyzed tube model of 0.5 m length with the cross-sectional area of 0.02 m × 0.02 m. It assigns frequency-dependent impedance BCs at the tube end. For FEM calculation, the tube model was discretized spatially by cubic FEs with edge length of 0.01 m. As an initial condition, the waveform of modulated Gaussian pulse with a peak frequency of 1 kHz was given to the vibration boundary at the tube inlet. The waveform and frequency characteristics of the used Gaussian pulse are presented in Figure 3. The sound pressures at two receiving points (x, y, z) = (0.44, 0.01, 0.01) and (0.45, 0.01, 0.01) were calculated using the 4th-E TD-FEM up to 1.0 s. Then the absorption coefficient and surface impedance of each absorber were computed using the transfer function method. As an accuracy measure for 4th-E TD-FEM, we defined the relative error of surface impedance $e_{\mathrm{Z}}$ between the theoretical values and the TD-FEM results as

$$e_{\mathrm{Z}}=\sqrt{\frac{1}{N_{\mathrm{f}}}\frac{{\sum }_{f=f_{\mathrm{l}}}^{f_{\mathrm{u}}}{|Z_{\mathrm{theory}}(f)-Z_{\mathrm{FEM}}(f)|}^2}{{\sum }_{f=f_{\mathrm{l}}}^{f_{\mathrm{u}}}{\left|Z_{\mathrm{theory}}(f)\right|}^2}}$$

(34)

where $N_{\mathrm{f}}$ is the total number of frequencies, and $f_{\mathrm{l}}$ and $f_{\mathrm{u}}$ are respectively the lowest frequency of 100 Hz and the highest frequency of 5 kHz. $Z_{\mathrm{theory}}(f)$ and $Z_{\mathrm{FEM}}(f)$ represent the surface impedance at frequency f calculated by the theory and the TD-FEM. $\Delta t$ was set to $\Delta t_{\mathrm{limit}}$. The convergence tolerance of the CR method was set to 10${}^{-6}$.

Figure 4a–c respectively present comparisons of normal-incidence absorption coefficient and normalized surface impedance among the theoretical results according to the transfer matrix method, the fitted results by the vector fitting method, and numerical results by the 4th-E TD-FEM for GW, NF, and MPP-GW. For the two porous-type sound absorbers of GW and NF, one can find good agreement at frequencies of 100 Hz to 5 kHz in both the absorption coefficient and surface impedance among the theory, fitted results, and 4th-E TD-FEM. The relative errors $e_{\mathrm{Z}}$s show respectively 0.035% and 0.011% for GW and NF. Those results show the validity of our formulation. We can also see good agreement in the absorption coefficient for the Helmholtz type sound absorber of MPP-GW. However, the impedance results showed that FEM and the fitted results have discrepancies from theory at frequencies higher than 3 kHz. Additionally, the larger relative error $e_{\mathrm{Z}}$ of 0.395% than those of GW and NF was observed. Therefore, at this stage, the present formulation still demands improvement for accurate modeling of MPP absorbers at higher frequencies. This is left as a subject for future work.

3.2. Small Cubic Room Problem

We examine further how accurate 4th-E TD-FEM is compared to FD-FEM. The latter is a standard selection for acoustics simulation because it can precisely reflect the sound absorption characteristics of sound absorbers. A fourth-order accurate FD-FEM [29] is used here to calculate reference solutions. Therefore, this examination can show FD-FEM users the potential of the proposed 4th-E TD-FEM for room acoustics simulation from the perspective of accuracy. Figure 5 shows the small, 1 m${}^3$ cubic room used for examinations. The room is a typical selection to verify the accuracy of wave-based room acoustics simulations. We assumed that the room has two sound-absorber configurations of Case 1 and Case 2, with their respectively different frequency-dependent impedance BCs settings. The sound absorber configuration of Case 1 has a GW ceiling absorber (red surface in Figure 5) with sound absorption characteristics shown in Figure 4a. The remaining boundaries have absorbing surfaces of normal-incidence sound absorption coefficient 0.052, which corresponds to normalized real-valued surface impedance of 74.4. Case 2 also has the GW ceiling absorber, and has NF floor absorber (blue surface in Figure 5) and MPP-GW absorber (green surface in Figure 5). The NF and MPP-GW absorbers have sound absorption characteristics shown in Figure 4b,c. The remaining boundaries are absorbing surfaces of normal-incidence sound absorption coefficient 0.052. A sound source was located at a corner (x, y, z) = (0, 0, 0). We used the same source signal as in Section 3.1 for 4th-E TD-FEM, whereas we used the volume acceleration of 1 ${\mathrm{m}}^3/{\mathrm{s}}^2$ for FD-FEM. Sound receivers were placed at six points on a line: $y=z=0.52$ with 0.2 m spacing. We calculated time responses up to 1 s at six receivers with a time interval $\Delta t_{\mathrm{limit}}$. Frequency responses were also calculated with a 1 Hz interval using FD-FEM. We used cubic FEs with a size of 0.04 m for spatial discretization.

The comparison was made using a mean sound pressure level (SPL) among six receivers at frequencies of 100 Hz to 1.2 kHz. Here, the SPL is a transfer function removing sound source characteristics in time-domain results according to a procedure described in an earlier report [45]. As a quantitative measure, the similarity of frequency responses between the FD-FEM result and the TD-FEM result were evaluated using the cross-correlate coefficient with respect to mean SPL $C{C}_{{\mathrm{L}}_{\mathrm{m}}}$ as

$$C{C}_{{\mathrm{L}}_{\mathrm{m}}}=\frac{{\sum }_{f=f_{\mathrm{l}}}^{f_{\mathrm{u}}}{L}_{\mathrm{m},\mathrm{FD}}(f)L_{\mathrm{m},\mathrm{TD}}(f)}{\sqrt{{\sum }_{f=f_{\mathrm{l}}}^{f_{\mathrm{u}}}{L}_{\mathrm{m},\mathrm{FD}}{(f)}^2}\sqrt{{\sum }_{f=f_{\mathrm{l}}}^{f_{\mathrm{u}}}{L}_{\mathrm{m},\mathrm{TD}}{(f)}^2}}$$

(35)

where $L_{\mathrm{m},\mathrm{FD}}(f)$ and $L_{\mathrm{m},\mathrm{TD}}(f)$ are respectively mean SPL calculated with the FD-FEM and TD-FEM at frequency f. Here, $f_{\mathrm{l}}$ and $f_{\mathrm{u}}$ were set to 100 Hz and 1.2 kHz, respectively. Hereafter, the tolerance value for the CR method was set to ${10}^{-4}$ which provided sufficient accuracy while maintaining efficiency. The effects of the tolerance value on the calculation results were examined in Appendix B.

Figure 6a,b respectively depict comparisons of mean SPLs between the FD-FEM and the 4th-E method with convergence tolerance of 10${}^{-4}$ for Case 1 and Case 2. In both cases, the TD-FEM results agree well with the FD-FEM results. The results revealed that the proposed 4th-E TD-FEM can conduct accurate acoustic simulation, including multiple frequency-dependent impedance BCs with a similar degree of accuracy as in FD-FEM in a qualitative sense. As a quantitative sense, the $C{C}_{{\mathrm{L}}_{\mathrm{m}}}$ values for Case 1 and Case 2 were respectively 0.998 and 0.999. In addition, we observed rapid convergence of CR method where numbers of mean iteration per time step were respectively 3.59 for Case 1 and 3.80 for Case 2.

4. Performance Comparison in Broadband Acoustic Modeling with Two Existing Implicit TD-FEMs

The present section presents a discussion of performances of the 4th-E TD-FEM and the Opt-E TD-FEM compared to the two existing implicit TD-FEMs [27,45], i.e., the fourth-order accurate implicit TD-FEM (4th-I) and the standard second-order implicit TD-FEM (2nd-I) in broadband room acoustic modeling. For Opt-E TD-FEM, we use the four optimized settings shown in Section 2.4. We will present a recommended setting for the Opt-E TD-FEM. The performance comparison was produced using the small cubic room problem with the sound absorber configuration of Case 1 used in Section 3.2. However, calculations were made at broad frequencies up to 2.5 kHz. All computations were performed using a workstation (Precision Tower 7810, Intel Xeon E5-2650 v4, 2.2 GHz; Dell Inc., Round Rock, TX, USA) with an Intel Fortran compiler (ver. 2019) and OpenMP parallel computation using eight threads.

4.1. Numerical Configuration

We assessed the performance achieved with three meshes M1, M2, and M3, each discretized with cubic FEs having different spatial resolutions. The coarsest mesh M1 has spatial resolution R of 5.5 elements per wavelength at 2.5 kHz, discretizing elements with edge size of 0.025 m. Mesh M2 has spatial resolution of $R=11$ using elements with 0.0125 m edge size. The finest mesh M3 has spatial resolution of $R=20.5$, employing elements with edge size of 0.0067 m. For performance examination, we use only the coarsest mesh M1 for 4th-E TD-FEM, Opt-E TD-FEM, and 4th-I TD-FEM because they are dispersion-reduced schemes that can accommodate larger elements than the standard 2nd-I TD-FEM. For standard 2nd-I TD-FEM, we use mesh M2, which satisfies a well-known rule of thumb for linear elements, i.e., at least ten elements per wavelength [67]. The finest M3 was used for computing the reference solution with 4th-I TD-FEM. The earlier report [45] describes that the 4th-I method can provide more accurate results than the fourth-order accurate FD-FEM with the same FE mesh in room acoustic modeling including the frequency-dependent impedance boundary. The presented scenario exhibits two aspects of the proposed methods: (1) how the 4th-E and Opt-E TD-FEMs perform well with the coarsest mesh compared to the standard 2nd-I TD-FEM, and (2) how the 4th-E and Opt-E TD-FEMs perform compared to 4th-I TD-FEM under the condition using the same coarsest mesh. Each method uses the stability-limit time interval $\Delta t_{\mathrm{limit}}$. Four optimization conditions are assumed for Opt-E TD-FEM, respectively minimizing the dispersion error at R = 6.87 (2 kHz), R = 6.25 (2.2 kHz), R = 5.5 (2.5 kHz), and R = 4.91 (2.8 kHz). The source point was set as (x, y, z) = (0, 0, 0). The six receivers were located on the line: $y=z=0.5$ with 0.2 m spacing. Each method calculated time responses up to 1 s using the source signal presented in Figure 3.

4.2. Measures of Accuracy

To quantify accuracy, we used the following measures. The first measure is cross-correlation coefficient $CC$, which evaluates similarity to the reference solution, defined as

$$CC=\frac{{\sum }_{f=f_{\mathrm{l}}}^{f_{\mathrm{u}}}{L}_{\mathrm{ref}}(f,{\mathit{r}}_i)L_{\mathrm{FEM}}(f,{\mathit{r}}_i)}{\sqrt{{\sum }_{f=f_{\mathrm{l}}}^{f_{\mathrm{u}}}{L}_{\mathrm{ref}}{(f,{\mathit{r}}_i)}^2}\sqrt{{\sum }_{f=f_{\mathrm{l}}}^{f_{\mathrm{u}}}{L}_{\mathrm{FEM}}{(f,{\mathit{r}}_i)}^2}}$$

(36)

Therein, $L_{\mathrm{ref}}(f,{\mathit{r}}_i)$ and $L_{\mathrm{FEM}}(f,{\mathit{r}}_i)$ are SPLs at frequency f on the receiver position of ${\mathit{r}}_i$ for the reference solution and TD-FEM results. They were calculated via discrete Fourier transformation to the resulting time responses. Here, the SPLs are transfer functions as used in the preceding section. Furthermore, $f_{\mathrm{l}}$ and $f_{\mathrm{u}}$ were respectively set to 100 Hz and 2.5 kHz. In an earlier study [45], $CC$ was used for time responses, but the present study applied $CC$ to frequency responses to evaluate similarity with respect to excitation with a flat spectrum.

We also defined accuracy measures using room acoustics parameters from practical aspects of architectural acoustics. Two room acoustics parameters were selected: early decay time ($EDT$) and $C_{50}$. They are, respectively, evaluated for perceived reverberance and clarity in speech. For the respective room acoustic parameters, we defined the following two error measures $e_{\mathrm{EDT}}$ and $e_{\mathrm{C}}$ as

$$e_{\mathrm{EDT}}(f_c)=\frac{1}{N_r}\sum _{i=1}^{N_r}\frac{|ED{T}_{\mathrm{ref}}(f_c,{\mathit{r}}_i)-ED{T}_{\mathrm{FEM}}(f_c,{\mathit{r}}_i)|}{ED{T}_{\mathrm{ref}}(f_c,{\mathit{r}}_i)}\times 100(\%)$$

(37)

$$e_{\mathrm{C}}(f_c)=\frac{1}{N_r}\sum _{i=1}^{N_r}|C_{\mathrm{ref}}(f_c,{\mathit{r}}_i)-C_{\mathrm{FEM}}(f_c,{\mathit{r}}_i)|\left(\mathrm{dB}\right)$$

(38)

Therein, $ED{T}_{\mathrm{ref}}(f_c,{\mathit{r}}_i)$ and $ED{T}_{\mathrm{FEM}}(f_c,{\mathit{r}}_i)$ respectively denote $EDT$ at 1/3 octave band center frequency $f_c$ at a receiver position ${\mathit{r}}_i$ for the reference solution and TD-FEM results. Also, $C_{\mathrm{ref}}(f_c,{\mathit{r}}_i)$ and $C_{\mathrm{FEM}}(f_c,{\mathit{r}}_i)$ represent those for $C_{50}$. Those parameters were calculated from the simulated impulse responses using ITA toolbox [68]. The reason for $EDT$ use instead of the reverberation time (RT) such as a $T_{20}$ is that the energy decay curve shows nonlinear decay with the unevenly distributed absorber placement [69].

4.3. Results and Discussion

First, we undertake comparisons of the SPL spectra at a receiver (x, y, z) = (0.6, 0.5, 0.5) among the reference solution, the Opt-E TD-FEM (optimized at R = 6.25), the 4th-E TD-FEM, 4th-I TD-FEM, and 2nd-I TD-FEM in Figure 7. The Opt-E method and 4th-I method show much better agreement with the reference solution in broadband than either the 4th-E method or the 2nd-I method. The 4th-E method shows discrepancies from the reference solution at frequencies higher than 2.1 kHz, where frequencies at which peaks and dips occur with shift to the lower frequency. The 2nd-I method using a finer mesh M2 shows the greatest discrepancies above 1.9 kHz and for frequencies at which peaks and dips shift to higher frequencies. These error characteristics of the 4th-E method and 2nd-I method are explainable by the theoretical dispersion error properties presented in Section 2.4, i.e., the dispersion error property of slower sound speed, as in the 4th-E method, engenders frequencies at which peaks and dips create a shift to a lower frequency than the reference solution. In comparison, the dispersion error property of faster sound speed, as in the 2nd-I method, produces the frequency shift to a higher frequency. Additionally, the much better accuracy of Opt-E and 4th-I TD-FEMs is explainable by their theoretical dispersion error properties presented in Figure 1, where these methods have theoretically much lower dispersion error than the 4th-E and 2nd-I TD-FEMs. Then, we show how the four optimizations work in the Opt-E TD-FEM. Figure 8 presents comparisons of SPLs above 1.8 kHz among the reference solution and the Opt-E TD-FEMs using the four optimizations. One can observe that each Opt-E TD-FEM shows excellent agreement around each optimized frequency: 2 kHz for R = 6.87, 2.2 kHz for R = 6.25, 2.5 kHz for R = 5.5, and 2.8 kHz for R = 4.91. In addition, at a frequency higher than the optimized frequency, frequencies at which peaks and dips occur with a shift to higher frequencies than the reference solution for all Opt-E methods. This error behavior can be explained by their theoretical dispersion error properties shown in Figure 1. One can suggest further that the optimization at too small R is ineffective, as might be apparent in the results for which R = 4.91.

Table 1. Comparison of $CC$ values among various TD-FEMs.

Method	4th-E	Opt-E (R = 6.874)	Opt-E (R = 6.25)	Opt-E (R = 5.5)
$CC$	0.843	0.938	0.945	0.918
Method	Opt-E ($\mathit{r}$ = 4.91)	4th-I	2nd-I
$CC$	0.835	0.964	0.802

shows the spatially averaged $CC$ values obtained using the respective methods. The 4th-I TD-FEM shows the highest CC value of 0.964. Three Opt-E TD-FEMs with R = 6.874, 6.25, and 5.5 are all high CC values larger than 0.91. As one might expect from the SPL comparison presented in Figure 7 above, Opt-E TD-FEM with R = 6.25 shows a comparable result to that obtained with 4th-I TD-FEM. Although Opt-E TD-FEM with R = 4.91 shows a slightly lower CC value than those for the other three Opt-E TD-FEMs, as expected from Figure 8, the value is still higher than the standard 2nd-I TD-FEM, which uses a resolution mesh that is two times higher. The 4th-E TD-FEM has a slightly higher CC value than the Opt-E TD-FEM with R = 4.91.

Table 2. Comparison of $e_{EDT}$ (%) among various TD-FEMs.

Method	Frequency, kHz
	1	1.25	1.6	2	2.5
4th-E	1.03	1.17	1.33	4.40	7.26
Opt-E (R = 6.874)	1.21	1.47	0.77	1.25	6.02
Opt-E (R = 6.25)	1.53	1.76	0.96	1.26	4.05
Opt-E (R = 5.5)	1.88	2.27	1.32	1.47	3.39
Opt-E (R = 4.91)	2.58	3.32	1.86	2.72	5.41
4th-I	1.30	1.49	1.10	2.36	9.98
2nd-I	3.15	4.63	1.88	3.72	10.26

and

Table 3. Comparison of $e_{\mathrm{C}}$ (dB) among various TD-FEMs.

Method	Frequency, kHz
	1	1.25	1.6	2	2.5
4th-E	0.04	0.06	0.29	0.65	1.66
Opt-E (R = 6.874)	0.10	0.19	0.26	0.39	1.26
Opt-E (R = 6.25)	0.12	0.22	0.26	0.40	1.00
Opt-E (R = 5.5)	0.16	0.26	0.25	0.54	0.47
Opt-E (R = 4.91)	0.19	0.32	0.29	0.67	0.76
4th-I	0.09	0.16	0.32	0.67	1.68
2nd-I	0.17	0.33	0.44	0.61	1.63

respectively present comparisons of $e_{\mathrm{EDT}}$ and $e_{\mathrm{C}}$ among the 4th-E TD-FEM, the Opt-E TD-FEMs with four optimized settings, 4th-I TD-FEM, and the 2nd-I TD-FEM. Herein, the results below 1 kHz were omitted because the error magnitude was smaller than the JND values for $EDT$ [70] and $C_{50}$ [71], i.e., 5% for $EDT$ and 1.1 dB for $C_{50}$. Only the two Opt-E TD-FEMs optimized at R = 6.25 and R = 5.5 produce smaller errors than JND values for both room acoustic parameters at all frequencies. It is particularly interesting that the two Opt-E TD-FEMs have higher accuracy at 2.5 kHz than the 4th-I TD-FEM. Considering the error estimations in the frequency and time domains, we recommend using the Opt-E TD-FEM optimized at R = 6.25 for room acoustic simulation at broadband.

As the end of this section, the computational costs among all TD-FEMs tested here are compared. It is noteworthy that the required memories for 4th-E and Opt-E TD-FEMs are equivalent, and that the computational complexities of both methods are also comparable because the mean iteration numbers of the CR method become almost identical to those shown in

Table 4. Mean numbers of iterations for explicit TD-FEMs.

Method	4th-E	Opt-E (R = 6.874)	Opt-E (R = 6.25)
Mean iteration	3.77	3.76	3.76
Method	Opt-E ($\mathit{r}$ = 5.5)	Opt-E ($\mathit{r}$ = 4.91)
Mean iteration	3.76	3.76

. Consequently, the computational costs of 4th-E and Opt-E TD-FEMs are presented as equivalent. The 4th-E and Opt-E TD-FEMs require 49.6 MBytes, which is slightly larger than 46.6 MBytes of 4th-I TD-FEM and which is 1/7.3 less than 359.0 MBytes of 2nd-I TD-FEMs. By contrast, 4th-E and Opt-E TD-FEMs present computational times of 82.1 s, which are respectively 1.4 times faster and 76.4 times faster than 116.9 s for the 4th-I method and 6272.2 s for the 2nd-I method. This computational cost comparison and the error estimation results described above suggest that the proposed Opt-E TD-FEM is a faster room acoustics solver than the existing two implicit TD-FEM solvers, while matching their accuracy.

5. Large-Scale Acoustics Simulation in an Auditorium Using a Highly Parallel Wave-Based Solver

This section presents a highly parallel Opt-E TD-FEM using the domain decomposition method for large-scale room acoustics simulation. Its applicability is demonstrated for an impulse response calculation in an auditorium model of 2271 ${\mathrm{m}}^3$, including approximate frequency components up to 3 kHz. Auditoriums are classical architectural spaces with large dimensions that need appropriate acoustical consideration for speech clarity. We selected this architectural space for the demonstration purpose of the proposed method because its acoustics simulation with the wave-based acoustics simulation methods is a challenging task due to the requirement of high computational cost. We use a supercomputer system consisting of 2000 nodes at Kyushu University: ITO, Subsystem A, Fujitsu Primergy CX2550/CX2560M4. Each node has two Intel Xeon Gold 6154 (3.0 GHz) with 18 cores. Also, Intel Fortran compiler ver. 2020 was used. The supercomputer has a distributed-memory type parallel computer among nodes, but shared-memory type parallel computation is available within each node.

5.1. Parallel Dissipation-Free and Dispersion-Optimized Explicit TD-FEM Using Domain Decomposition

The domain decomposition method (DDM) is a general term for division of the computational domain into multiple subdomains and performance of calculations in the respective subdomains to obtain a solution. In the parallel computation method based on DDM, each node or processor is assigned to compute each subdomain and to execute the computation. Using the DDM-based parallel computation technique, all computational processes in the Opt-E TD-FEM, i.e., matrices construction processes and time-marching scheme, can be parallelized. Communications using MPI are necessary when using data stored for each node.

For DDM-based parallel computations, one must first partition an FE mesh to subdomain meshes that include communication information with a partitioning algorithm. This study uses Metis Library [72] for the partitioning process, and the FE mesh is divided into non-overlapping subdomain meshes. With subdomain meshes, the finite element matrices construction processes, i.e., constructions of $\mathit{M}$, $\mathit{K}$, and $\mathit{D}$, are performed in parallel at each node or processor. For ease of implementation, the surface matrix $\mathit{C}$ is assembled for the present study with the nonpartitioned surface mesh. Then, it is sorted into each process. After finishing the matrix construction processes, the time-marching computation consisting of Equations (22), (26) and (33) is also performed in parallel at each node or processor with communications using MPI. Communications appear in the following processes: (1) the assembly process of $\mathit{D}$, which performs only once before the time-marching computation, (2) the two sparse matrix–vector product operations per time step of $\mathit{Mv}$ in Equation (22) and $\mathit{Kp}$ in Equation (33), and the solution process of Equation (26) using the CR method, which includes sparse matrix–vector operations and inner product operations. Thread-parallelization is available within a shared-memory type architecture. We also use thread-parallelization using OpenMP.

As a detailed calculation procedure, the sparse matrix-vector product operations $\mathit{Mv}$ can be described as the summation of sparse matrix-vector product results at each subdomain.

$$\mathit{Mv}=\sum _{i=1}^{N_{\mathrm{sd}}}{\mathit{M}}_i{\mathit{v}}_i=\sum _{i=1}^{N_{\mathrm{sd}}}{\mathit{w}}_i$$

(39)

where $N_{\mathrm{sd}}$ is the total number of subdomains, ${\mathit{M}}_i$ and ${\mathit{v}}_i$ are respectively $\mathit{M}$ and $\mathit{v}$ in the i-th subdomain, and ${\mathit{w}}_i$ is the resulting vector from the matrix-vector product in the i-th subdomain. In the DDM-based parallel computation, the sparse matrix-vector product operations are first performed at each subdomain and the vector ${\mathit{w}}_i$ of each subdomain are computed as a result. Then, the vector ${\mathit{w}}_i$ of each subdomain has to be added the components of nodes at overlapping boundary interfaces of other subdomains, with a communication. This treatment is also necessary for constructing the matrix $\mathit{D}$. Furthermore, one must consider the influence of duplicate boundary interfaces among subdomains to calculate the inner product, i.e., the vector component at the boundary interface’s nodes must be divided by the number of subdomains to which the nodes belong. With this consideration, the inner product is computed as

$$(\mathit{q}·\mathit{q})=\sum _{i=1}^{N_{\mathrm{sd}}}({\mathit{q}}_i·{\mathit{q}}_i/{\mathit{N}}_{\mathrm{bd}})=\sum _{i=1}^{N_{\mathrm{sd}}}{s}_i$$

(40)

where $\mathit{q}$ is arbitrary vectors, ${\mathit{q}}_i$ is $\mathit{q}$ in the i-th subdomain, ${\mathit{N}}_{\mathrm{bd}}$ is the vector stored the number of subdomains to which the nodes belong, and $s_i$ is the resulting scalar from the inner product in the i-th subdomain. The inner product is first computed at each subdomain, and then the resulting scalar $s_i$ of each subdomain is summed via a communication.

5.2. Room Configuration and Boundary Conditions

Figure 9 presents the analyzed auditorium model. We analyzed an auditorium with two BCs (Cond. 1 and Cond. 2). In Cond. 1, the red, blue, and green surfaces in Figure 9 were, respectively, assigned the frequency-dependent surface impedance of GW, NF, and MPP-GW used in Section 3. Other surfaces were reflective boundaries with real-valued frequency-independent impedance corresponding to the normal-incidence absorption coefficient of 0.1. In Cond. 2, 12 additional sound-absorbing panels with GW absorption characteristics were further installed into the yellow surfaces on sidewalls as sidewall absorbers to assess effects of sidewall absorbers in a large auditorium. The sound-absorbing panels have size of 1.74 m × 1.2 m. It is noteworthy that the used BCs are a virtual scenario to demonstrate the applicability of parallel Opt-E TD-FEM to practical applications. The source point was placed on (x, y, z) = (2.1, 0, 1.2). We set 14 sound receivers at the coordinates presented in

Table 5. Locations of receivers in the auditorium model.

Receiver	(x, y, z)	Receiver	(x, y, z)
R1	(8.07, 0, 0.9)	R8	(13.3, 4.6, 2.7)
R2	(8.07, 4.6, 0.9)	R9	(15.0, 0, 3.3)
R3	(9.81, 0, 1.5)	R10	(15.0, 4.6, 3.3)
R4	(9.81, 4.6, 1.5)	R11	(16.8, 0, 3.9)
R5	(11.6, 0, 2.1)	R12	(16.8, 4.6, 3.9)
R6	(11.6, 4.6, 2.1)	R13	(18.5, 0, 4.5)
R7	(13.3, 0, 2.7)	R14	(18.5, 4.6, 4.5)

.

5.3. Simulation Outline

For acoustics simulations, the auditorium model was spatially discretized using variously shaped 8-node hexahedral FEs having a maximum edge size smaller than 0.025 m. The resulting FE mesh has 146,583,255 DOFs. For DDM-based parallel computation, we partitioned the FE mesh to 256 subdomain meshes using Metis Library [72]. Figure 10 presents the image of decomposed auditorium model with 256 subdomains. Also, Figure 11 shows a workload distribution among 256 subdomains, i.e., DOF of each subdomain and number of communication data per communication in each subdomain. The average values for DOF and communication data among all subdomains were respectively 596,287 and 48,451. The coefficient of variance for the distributed DOF is 0.72%, meaning that the DOF of the global auditorium model was evenly arranged to each subdomain. In addition, the coefficient of variance for communication data among all subdomains is 18.2%, which was larger than that for DOF. However, a workload-balanced domain decomposition among the subdomains without the extreme deviation of the partitioning was achieved. Because we use two thread-parallelization for each node or processor, parallel computation using 512 cores was performed, i.e., 256 process-parallelization by MPI and two thread-parallelization by OpenMP. For effective dispersion optimization with FE mesh composed of nonuniform elements, each element sets the integration points to minimize dispersion error at 2.2 kHz with the maximum edge length per element. Actually, this optimization procedure is based on the earlier work [37]. Considering the stability condition of Equation (32) using the smallest edge size in the FE mesh, a time interval was set as 1/30,000 s in both BCs. We calculated impulse responses up to 3 s with the Gaussian pulse sound source of Figure 3. Using the computed impulse responses, we made comparisons of the impulse responses (IRs), energy decay curves (EDCs), and four room acoustic parameters (reverberation time $T_{30}$, $EDT$, $C_{50}$ and sound strength G) between the two conditions. The EDCs and room acoustic parameters were estimated according to ISO3382-1 [70]. The room acoustic parameters were calculated for the 125 Hz to 2 kHz octave band. In addition, the parallel performance of the proposed solver was examined by evaluating the speedup against the non-parallelized case. The evaluation is made with analyzed time up to 0.1 s (3000 steps) because the computational time of the non-parallelized solver exceeds the resource limit on the supercomputer when calculating 3 s impulse response.

5.4. Results and Discussion

First, we can show the excellent computational speed of the proposed parallel Opt-E TD-FEM. The computational times for IR calculations up to 3 s were only 8878 s and 8996 s, respectively, for Cond. 1 and Cond. 2. The slight increase of computational time for Cond. 2 occurs because the condition includes more frequency-dependent impedance BCs for the GW sound-absorbing panels. The speedups against the non-parallelized solver are 336 and 341 for Cond. 1 and Cond. 2, respectively, showing its notable parallel performance. This result demonstrates that the practicality of our proposed method in room acoustics design. It is noteworthy that the computational speed can be enhanced further using more CPU cores. At this stage, our computational code is only insufficiently optimized. A detailed performance study of massively parallel computing environments remains as a subject for future work.

Figure 12a–d, respectively present the simulated sound field on the $zx$ plane (y = 0) at 5 ms, 12.5 ms, 20 ms, and 27.5 ms. In fact, the figure shows only the sound field of Cond. 1. There is no difference between the two cases at the snapped times because the sound does not reach the sidewalls, where different BCs were assumed. Isotopic sound wave propagation can be observed, indicating that the proposed Opt-E TD-FEM reduces the dispersion error effect. Figure 12c,d present diffracted sound around the seating area. Additionally, the sound absorption effect of the ceiling absorber is evident in Figure 12d, where the reflected sound from the ceiling (green wavefront) was weakened more than other sounds (yellow wavefront). This visualization presents an intuitive understanding of the sound absorber effectiveness.

Figure 13a,b portray IRs and EDCs at a receiver R13 for Cond. 1 and Cond. 2. An acoustic defect called flutter echo is apparent for the Cond. 1 result. The reason for these multiple reflections is the low absorptivity of sidewalls. The time interval of 44.2 ms between strong reflections corresponds to the sound travel time between the side walls. One can also find multiple reflections effects in the EDCs, where nonlinear decay is observed at all frequencies. Especially, the EDCs at 1 kHz and 2 kHz are curved to a considerable degree because a high sound absorptive ceiling absorber gives only a small contribution to sound waves with a shorter wavelength propagating parallel to the absorber surface. The multiple reflections can be alleviated by the additional sound-absorbing panels as in Cond. 2 results. Furthermore, strong reflected sounds in early times were reduced. One can see smoother EDCs at frequencies higher than 500 Hz than those in Cond. 1. The decay properties at 1 kHz and 2 kHz were improved much more than at 500 Hz, which reflects the absorption property effects of additional sound-absorbing panels.

The acoustics improvement effects achieved using sound-absorbing panels can be observed further in four room acoustical parameters. Figure 14a–d, respectively show spatially averaged values of (a) $T_{30}$, (b) $EDT$, (c) $C_{50}$, and (d) G for Cond. 1 and Cond. 2. The figures also present the standard deviation over all receivers’ results in each acoustic parameter as an error bar. Figure 14a shows that $T_{30}$ in Cond. 1 decreases up to 500 Hz, and then it increases with frequency. This increase results from the above-explained effect of the nonlinear decay in EDCs. However, $T_{30}$ is reduced by the additional sound-absorbing panel. One can find $T_{30}$ in Cond. 2 decreases to about 1.3 s at 1 kHz and 2 kHz. Additionally, the spatial variances of $T_{30}$ at high frequencies were reduced thanks to the additional absorbers. The coefficient of variance in Cond. 2 is reduced from 12.6% to 6.4% at 1 kHz compared to the results in Cond. 1. At 2 kHz, it is reduced from 8.6% to 3.0%. That is because the effect of flutter echo differs on the receiver positions in Cond. 1. The other room acoustical parameters $EDT$ and G at 1 kHz and 2 kHz and $C_{50}$ at 2 kHz in Cond. 2 are also changed considerably from Cond. 1 by more than each JND value, i.e., 5% for EDT, 1.1 dB for $C_{50}$ and 1 dB for G [70]. Those changes become greater with frequency. The results presented in this section demonstrate clearly that practicality of the proposed parallel Opt-E TD-FEM for practical room acoustic simulations including multiple sound absorbers.

6. Conclusions

As the most notable contribution of this study, we proposed a DDM-based parallel wave-based acoustics simulation method using dissipation-free and dispersion-optimized explicit TD-FEM for large-scale room acoustics problems. The proposed parallel solver can accommodate the sound absorption effects of various sound absorbers as the locally reacting frequency-dependent impedance BCs. This study used three numerical problems to assess its basic performance: an impedance tube problem, a small cubic room problem, and a large auditorium problem. The assessment includes (1) accuracy comparison with FD-FEM, which is a standard selection for acoustics simulation, and (2) performance comparison with existing TD-FEMs. It is noteworthy that the dissipation-free and dispersion-optimized explicit TD-FEMs abbreviated to the 4th-E and Opt-E TD-FEMs themselves have been proposed in the earlier work [37] for acoustics simulation with rigid boundaries. The novelty of the present study against earlier work is the proposal of a time-marching scheme that is able to address multiple frequency-dependent impedance BCs and the proposal of a DDM-based parallel solver for large-scale problems.

We first validated the proposed method using the impedance tube problem and a small cubic room problem, particularly addressing the 4th-E TD-FEM. The impedance tube problems with three sound absorbers of GW, NF, and MPP-GW demonstrated overall that the 4th-E TD-FEM can simulate their sound absorption characteristics accurately compared to theoretical results obtained using transfer matrix method. However, we found a future subject presenting difficulty: reproduction of the complex-valued surface impedance of MPP-GW at high frequencies. Small cubic room problems with multiple frequency-dependent impedance BCs that examines the accuracy of the 4th-E TD-FEM compared to the fourth-order accurate FD-FEM revealed that the 4th-E TD-FEM has a similar level of accuracy as that of FD-FEM.

Subsequently, we investigated the performance gain of the 4th-E TD-FEM and the Opt-E TD-FEM with four optimized setups against two existing implicit TD-FEMs: 4th-I TD-FEM and 2nd-I TD-FEM. Results demonstrated the outstanding computational speed of the Opt-E TD-FEM using the optimization at $R=6.25$, which is 1.4 times faster than the 4th-I TD-FEM and 76.4 times faster than the 2nd-I TD-FEM while maintaining equivalent accuracy.

Finally, we tested the parallel computational performance of the DDM-based parallel Opt-E TD-FEM using a large auditorium model of 2271 ${\mathrm{m}}^3$ with multiple sound absorbers. The auditorium model has about one hundred fifty million DOFs. The results revealed that the parallel Opt-E TD-FEM can simulate impulse responses with a time length of 3 s, including frequency component up to 3 kHz within 9000 s under the use of 512 CPU-cores with MPI-OpenMP hybrid parallel computational framework. In addition, the parallel solver presented the notable speedup against the non-parallelized solver, with 336–341 times faster computational speed. To conclude, the present parallel wave-based solver has high potential for use with practical room acoustics design tools. It is also helpful for FD-FEM users because this method can be coded easily from standard linear FD-FEM code with a simple modification of integration points for element matrices and implementation of the three-step time integration method.