pffdtd: wave-based FDTD for room acoustics

This is the technical note of the pffdtd fork ST-LINE maintains: what a wave-based FDTD solver does for room acoustics, what it costs, where it goes wrong, and what the fork changed with respect to Brian Hamilton’s original — the numerical scheme no, the CUDA engine yes. The repository is at github.com/stefanofante/pffdtd; the showcase card is in Open Lab.

FDTD for room acoustics

The Finite-Difference Time-Domain (FDTD) method solves linear acoustics in the time domain by discretising directly the first-order equations that link pressure and particle velocity. With no mean flow and for small perturbations, conservation of momentum and of mass give the system

$$\rho \frac{\partial \mathbf{v}}{\partial t}=-\nabla p,\frac{1}{\rho c^2}\frac{\partial p}{\partial t}=-\nabla \cdot \mathbf{v}$$

where $p$ is the acoustic pressure, $\mathbf{v}$ the particle velocity, $\rho$ the air density and $c$ the speed of sound. Eliminating the velocity yields the scalar wave equation ${\partial }_{tt}p=c^2{\nabla }^{2}p$, and it is this that FDTD marches step by step on a regular grid, replacing the derivatives with centred finite differences in time and space.

Cartesian vs FCC stencil

The choice of stencil — the set of neighbouring nodes that approximate the Laplacian — governs accuracy and cost. The simplest scheme is the 7-point stencil on a Cartesian grid (the central node plus its six nearest axial neighbours). It is cheap but strongly anisotropic: the numerical wave speed depends on the propagation direction relative to the grid, and this numerical dispersion introduces phase error that grows with frequency.

A face-centred cubic (FCC) grid with a 13-point stencil distributes the neighbours more isotropically. For a given tolerated dispersion, FCC allows a coarser spatial step and ends up costing about five times less memory than the 7-point Cartesian stencil — a decisive factor, because in FDTD the node count grows with the cube of the maximum frequency and RAM consumption is the dominant constraint.

Stability, CFL and energy conservation

An explicit time-marching scheme is stable only if the time step respects the Courant–Friedrichs–Lewy (CFL) condition: information may not cross more than one cell per step. Compactly,

$$\lambda =\frac{c\Delta t}{\Delta x}\le {\lambda }_{\max }$$

with ${\lambda }_{\max }$ depending on the stencil (it is $1/\sqrt{3}$ for the 3-D 7-point Cartesian scheme). Exceeding the limit blows the solution up; staying just below it maximises per-node accuracy.

A deeper property is energy conservation. The scheme can be cast in passive form, with a discrete energy functional that never grows in the absence of sources and dissipation. In double precision this conservation is verifiable to machine precision: the energy balance closes to within the numerical epsilon step after step, and that is the most stringent sanity check a wave-based solver can offer. In single precision conservation holds only within a wider margin, so the engine adopts safeguards (energy monitoring, careful operation ordering) to avoid slow drift in long runs.

Frequency-dependent impedance boundaries (ADE)

Real materials do not have constant impedance: they absorb differently at different frequencies. Modelling a boundary as a frequency-dependent impedance condition $Z(\omega )$ means, in the time domain, a convolution between pressure and the boundary response — costly and long-memoried. The auxiliary differential equation (ADE) technique sidesteps the convolution: it represents the impedance as a rational function of $\omega$ and introduces, for each boundary node, a few auxiliary variables that evolve according to a local ODE. The global convolution thus becomes a small system of ODEs updated in place, compatible with the explicit, local nature of FDTD.

What FDTD costs, in numbers

That the node count grows with the cube of frequency and that RAM is the dominant constraint is true, but while it stays qualitative it does not say which problems are tractable. Put it in figures on a 500 m³ room, with six points per wavelength at the maximum frequency, 1.5 s of simulated tail, a 7-point Cartesian stencil (CFL 1/√3) and three fields in double precision, that is 24 bytes per node:

The law is simple and brutal. The spatial step scales as 1/f, so nodes as f³; but CFL ties the time step to the spatial one, so time steps also grow as f. Total work therefore goes as f⁴: going up one octave costs eight times the memory and sixteen times the time. Between 250 Hz and 8 kHz, five octaves, there are six orders of magnitude of memory and six of work.

This is what makes FDTD a method for the low end of the spectrum, not a stylistic choice.

What the FCC stencil buys

Put that way, the FCC’s factor of five in memory looks like one optimisation among many. In fact, inside a cubic law, five times the memory is worth 5^(1/3) = 1.71 times the maximum frequency, that is about three quarters of an octave more bandwidth on the same machine. And since a spatial step 1.71 times coarser also brings 1.71 times fewer time steps, the saving in total work is 8.5-fold, not fivefold.

Three quarters of an octave is the difference between stopping at 5 kHz and reaching 8.5 kHz, or between a run that fits in 32 GB and one that asks for 160. Put in those terms it is clear why the choice of stencil comes before any micro-optimisation of the kernel.

Where the ceiling sits, in practice

The bandwidth ceiling depends on volume, and there too the dependence is cubic. Under the same assumptions, the maximum frequency that fits in memory:

This should be compared with the frequency below which a wave method is genuinely necessary, that is the Schroeder frequency, where the room modes separate and the statistical description fails (see the reverberation time wiki): 179 Hz for the 50 m³ control room at T60 0.4 s, 89 Hz for the hall at T60 1 s, 23 Hz for the concert hall at T60 2 s.

The comparison is instructive, and runs opposite to what one would expect: the FDTD ceiling sits one or two decades above the frequency at which a wave method becomes indispensable. On an ordinary machine, wave methods cover far more spectrum than diffuse-field physics requires.

Hence the real reason for hybrid stacks, which is not an inability to reach f_S. It is that (a) you want wave accuracy well above f_S, where early reflections and diffraction matter and geometrical acoustics is crude precisely there; and (b) in inverse design the forward simulation is not run once but hundreds of times. A 2 kHz run on the 500 m³ room costs 6.7 · 10¹¹ node updates: an optimisation with a hundred iterations and two solves per gradient asks for 1.3 · 10¹⁴. It is in that regime that cost per degree of freedom becomes the parameter deciding everything — and it is why the next section arrives at discontinuous Galerkin.

Staircasing and its correction

A regular grid represents axis-aligned walls well, but an inclined or curved surface is approximated as a staircase (staircasing). The problem is not merely cosmetic: the effective boundary area seen by the scheme is systematically wrong — the steps inflate the exposed area — and since total absorption is proportional to area, staircasing leads to overestimating absorption and hence, at times, to mis-estimating decay times. On unfavourable geometries the reverberation-time error can reach around 50 %.

The effective-area correction fixes this by weighting each boundary node’s contribution with the area actually intercepted by the continuous surface, rather than with the staircase area. Applied on progressively finer grids, it brings the error below 1 %. It remains a remedy, though: the underlying geometry is still staircased, and this is one of the structural reasons why, when exact geometry is needed, one moves to a body-conforming method (see Discontinuous Galerkin below).

The pffdtd fork: what was optimised

pffdtd originates as the FDTD simulator by Brian Hamilton (University of Edinburgh, 2021), released under the MIT licence. Our fork (github.com/stefanofante/pffdtd) leaves the numerical scheme unchanged: for identical inputs, the fork’s forward output matches the reference Python engine to machine precision. We did not touch the physics; we worked underneath it, on the CUDA GPU engine, where the original implementation left performance and portability on the table.

The optimisations, in brief:

• Native-architecture builds. Instead of compiling for sm_35 (Kepler) and relying on runtime PTX-JIT recompilation on modern GPUs, the fork builds for the card’s actual architecture. This removes the first-launch JIT latency and the inefficiencies of code generated for an obsolete ISA.
• Real multi-GPU peer access. On problems exceeding a single card, the original silently staged data through host RAM, bottlenecking on the PCIe bus. The fork enables direct peer access between GPUs, when topology allows, so domain partitions communicate device-to-device.
• Batched source injection. The original launched a <<<1,1>>> kernel — a single thread — per source and per timestep, wasting the GPU’s parallelism. The fork injects sources in batch, in a single kernel.
• Block-wise read-out. Receiver extraction happens in blocks rather than point-by-point, reducing the number of transfers.
• Runtime-queried memory budget. Rather than a hardcoded threshold, the engine asks the GPU how much memory is available (cudaMemGetInfo) and sizes the problem accordingly — it discovers the hardware and adapts, instead of assuming.

Two hardware tuning profiles

Tuning was done on two deliberately different machines: an RTX 4500 Ada with 24 GB of dedicated VRAM and a DGX Spark GB10 with 128 GB of unified memory. The difference in memory behaviour is instructive and drove the sizing logic. On dedicated VRAM, over-allocation fails hard: once capacity is exceeded the allocation errors out and the run stops. On unified memory the system degrades gracefully, paging between GPU and host: the run continues, but slows down. Querying the budget at runtime and sizing the domain below the real threshold avoids both pathologies — the hard crash on one side, the silent degradation on the other.

References

• B. Hamilton, pffdtd — FDTD solver for room acoustics, 2021, MIT licence. github.com/bsxfun/pffdtd.
• F. Mondet et al., few-parameter fractional impedance model for acoustic surfaces, 2020. DOI 10.1016/j.apacoust.2019.04.034.
• J. S. Hesthaven, T. Warburton, Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications, Springer. DOI 10.1007/978-0-387-72067-8.
• C. F. Eyring, Reverberation time in “dead” rooms, J. Acoust. Soc. Am., 1930. DOI 10.1121/1.1915175.
• M. R. Schroeder, on the transition frequency between the modal and diffuse regimes. DOI 10.1121/1.1909343.
• L. Aspöck et al., BRAS — Benchmark for Room Acoustical Simulation, TU Berlin, 2020. DOI 10.14279/depositonce-6726.3.

The reported values come from the project documentation and the fork’s README; DOIs not directly verified are cited by author, year and title without a DOI.