## 4 Finite-Difference Time-Domain Simulations  
In room acoustics, multiple simulation methods aim to directly solve the wave equation in homogeneous form:  
\[\left(\
abla^{2}-\frac{1}{c_{0}^{2}}\frac{\partial^{2}}{\partial t^{2}}\right)p=0 \tag{17}\]  
Figure 12: Sound absorption coefficient and normalized surface impedance for four DL configurations. (**a**) Configuration A. (**b**) Configuration B. (**c**) Configuration C. (**d**) Configuration D.  
where \(\
abla^{2}\) is the 3D Laplacian operator, \(c_{0}\) is the speed of sound at \(T=20\,^{\circ}\)C and RH\(=50\%\), and \(p=p(\mathbf{x},t)\) is the pressure (deviation from ambient). The finite-difference time-domain (FDTD) method is among the oldest numerical methods to solve time-dependent partial differential equations (PDEs), like the wave equation [36]. In FDTD methods, space and time continuous domains are typically discretized with regular Cartesian grids and updating recursions occur in the nodal points to calculate the acoustic properties such as the pressure \(p(\mathbf{x},t)\) or the particle velocity potential \(u(\mathbf{x},t)\)[37]. The solution \(p(\mathbf{x},t)\) of the wave equation, with \(\mathbf{x}\in\mathbb{R}^{3}\), is approximated by a grid function \(p_{l,m,p}^{u}\) at spatiotemporal points \(x=lh,y=mh\), \(z=ph\) and \(t=nk\), with \(l\), \(m\), \(n\) and \(p\) representing integer numbers, \(h\) the grid spacing, and \(k\) the time step [38, 39]. A large variety of explicit FDTD methods follows the same general scheme:  
\[\delta_{t}^{2}p_{l,m,p}^{n}=\lambda^{2}[(\delta_{x}^{2}+\delta_{y}^{2}+\delta_ {z}^{2})+a(\delta_{x}^{2}\delta_{y}^{2}+\delta_{x}^{2}\delta_{z}^{2}+\delta_{ y}^{2}\delta_{z}^{2})+b(\delta_{x}^{2}\delta_{y}^{2}\delta_{z}^{2})]p_{l,m,p}^{n} \tag{18}\]  
where \(\lambda\) is the dimensionless quantity defined as the Courant number \(\lambda=ck/h\) and \(a\) and \(b\) are the specific coefficients of a general family of compact explicit schemes [40]. The operators \(\delta_{t}^{2}\) and \(\delta_{x}^{2}\) act on \(p_{l,m,p}^{n}\) as follows:  
\[\delta_{t}^{2}p_{l,m,p}^{n}=p_{l,m,p}^{n+1}-2p_{l,m,p}^{n}+p_{l,m,p}^{n-1}, \qquad\delta_{x}^{2}p_{l,m,p}^{n}=p_{l+1,m,p}^{n}-2p_{l,m,p}^{n}+p_{l-1,m,p}^{ n} \tag{19}\]  
and similarly for \(\delta_{y}^{2}\) and \(\delta_{z}^{2}\). The value of the Courant number is closely correlated to the stability condition of the system, which is guaranteed avoiding the exponential growth of the numerical solution (see, e.g., [41]). In the simple case of the seven-point scheme (a = b = 0), the stability condition is expressed as:  
\[\lambda\leq\frac{1}{\sqrt{3}} \tag{20}\]  
which is the so-called Courant-Friedrichs-Lewy (CLF) condition [36].  
However, in this paper the cubic close-packed (CCP) scheme (see Figure 13), with \(a=0.25\), \(b=0\) and \(\lambda=1\), is employed for its favorable numerical dispersion properties [41]. The consequent maximal time step, for a certain grid spacing, is k < h/c. As is the case for any 3D FDTD scheme, the computational cost is proportional to \(h^{-3}\), which can quickly lead to large simulation grids. Then, such grid-based simulations are natural candidates for parallel acceleration on graphics processing units (GPUs) [42].