## 2 Theory
### Spatial Discretization with Finite Elements  
The inhomogeneous lossy Helmholtz equation of Equation (4) is solved numerically with standard Galerkin FEM discretization to compute the complex-valued sound pressure, including the air absorption effect. The weak form of Equation (4) is  
\[\int_{\Omega}\
abla\chi(\
abla\hat{p})^{\text{T}}\text{d}\Omega-\int_{\Gamma} \chi\frac{\partial\hat{p}}{\partial n}\text{d}\Gamma-\vec{k}^{2}\int_{\Omega} \chi\hat{p}\text{d}\Omega=\text{j}\omega\rho_{0}\hat{q}\int_{\Omega}\chi\hat{ \delta}(\mathbf{r}-\mathbf{r}_{\text{s}})\text{d}\Omega, \tag{6}\]  
where \(\Omega\) denotes the domain, \(\Gamma\) stands for the boundary of \(\Omega\), and \(\chi\) is the weight function. In FEM, sound pressure \(\hat{p}(\mathbf{r})\) at point \(\mathbf{r}\) within element \(\Omega_{e}\) is approximated by a linear combination of two vectors \(\mathbf{N}\) and \(\mathbf{p}_{e}\) as \(\hat{p}(\mathbf{r})=\mathbf{N}\mathbf{p}_{e}\), where \(\mathbf{N}\) and \(\mathbf{p}_{e}\), respectively, denote the shape function vector and nodal sound pressure vector. The following Lagrange polynomials are used for \(\mathbf{N}\):  
\[l_{a}(\xi^{m})=\prod_{b=1,(b\
eq a)}^{p^{\xi^{m}}+1}\frac{\xi^{m}-\xi^{m}_{b}} {\xi^{m}_{a}-\xi^{m}_{b}}, \tag{7}\]  
where \(l_{a}(\xi^{m})\) is the Lagrange polynomial for the \(\xi^{m}\) direction in the local coordinate system. Subscript \(a\) denotes the node location, \(\xi^{m}_{a}\) is the corresponding local coordinate, \(p^{\xi^{m}}+1\) is the node number used for FE approximation, and \(p^{\xi^{m}}\) is the order of the polynomial. This study specifically addresses 2D room acoustics simulation discretized using four-node linear quadrilateral elements (Q4s) with \(p^{\xi^{1}}=p^{\xi^{2}}=1\). In this case, the shape function vector \(N_{i}(i=1,2,3,4)\) components are given as  
\[N_{i}(\xi^{1},\xi^{2})=l_{a}(\xi^{1})l_{a}(\xi^{2}). \tag{8}\]  
With Q4s, \(\mathbf{N}\) and \(\mathbf{p}_{e}\) are given, respectively, as \(\mathbf{N}=[N_{1},N_{2},N_{3},N_{4}]\) and \(\mathbf{p}_{e}=[\hat{p_{1}},\hat{p_{2}},\hat{p_{3}},\hat{p_{4}}]^{\text{T}}\). The same function form as \(\hat{p}(\mathbf{r})\) is applied to the weight function for the spatial discretization of Equation (6). Using sound field approximation within \(\Omega_{e}\) and inserting three boundary conditions (rigid boundary \(\Gamma_{r}\), vibration boundary \(\Gamma_{v}\), and impedance boundary \(\Gamma_{a}\) conditions) to Equation (6) leads to the following linear system of equations as  
\[(\mathbf{K}-\vec{k}^{2}\mathbf{M}+\text{j}\omega\mathbf{C})\mathbf{p}=\mathbf{f}, \tag{9}\]  
where \(\mathbf{K}\), \(\mathbf{M}\), and \(\mathbf{C}\), respectively, represent the stiffness, mass, and dissipation matrices, and where \(\mathbf{f}\) is the external force vector. Sound pressure vector \(\mathbf{p}\) composed of complex-valued sound pressure at all nodes is computed by solving Equation (9) using a linear system solver, such as a direct solver [50] or iterative solver [51]. Since the coefficient matrix of Equation (9) becomes a complex-valued sparse symmetric matrix with many zero components, the use of sparse linear system solvers is helpful to conserve memory when predicting a large-scale sound field with vastly numerous degrees of freedom (DOF). Thisreport uses two sparse linear system solvers. For the accuracy assessment in Section 3, a sparse direct solver called PARDISO (included in the Intel Math Kernel Library) is used because sparse direct solvers produce accurate solutions. Additionally, an iterative solver called the CSQMOR method [52] with a diagonal scaling preconditioning is used as a linear system solver for the office problem in Section 4 for efficiency estimations because iterative solvers can have the potential to solve large-scale problems with lower memory usage and with less computational time than direct solvers. The convergence tolerance is set as \(10^{-4}\).