Simulating 3D Room Modes: Voxel Finite-Difference Helmholtz Solvers for the Browser

TL;DR: Ray tracing breaks down in small room acoustics below 200 Hz because sound behaves as discrete standing waves rather than geometric rays. We built roomtreatment.diy by discretizing arbitrary 3D room geometries into a uniform voxel grid ($h = 0.15\text{ m}$ pitch, roughly $\lambda/11$ at 200 Hz, coarsened automatically for large rooms), solving the spatial negative Laplacian $-\nabla^2$ on the server with a SLEPc Krylov–Schur eigensolver (MUMPS shift-invert), and applying first-order acoustic perturbation theory to predict modal RT60 decay and speaker/listener nulls per band.

Most web-based audio tools treat room acoustics like optics. They draw straight lines from a speaker to a wall, bounce them off at equal angles, and call it a day.

That works fine for high frequencies. But below 200 Hz, ray tracing is a lie.

In a typical home studio or listening room, a 50 Hz sound wave is nearly 7 meters long. It doesn't bounce around like a laser beam. It fills the room, wraps around itself, and creates massive acoustic standing waves—room modes. You get localized peaks where bass booms uncontrollably and deep nulls where low frequencies vanish entirely.

If you want to simulate how bass actually behaves in real, non-shoebox rooms (L-shapes, slanted ceilings, alcoves), you have to solve the 3D wave equation directly. Here is how we implemented a server-side Finite-Difference Helmholtz eigensolver that feeds an interactive WebGL viewer.

The Physical Model: Helmholtz Equation & Boundary Conditions

In the low-frequency domain ($f < 200\text{ Hz}$), acoustic sound pressure $p(\mathbf{r}, \omega)$ inside a room volume $\Omega$ satisfies the inhomogeneous Helmholtz equation:

$\nabla^2 p(\mathbf{r}) + k^2 p(\mathbf{r}) = -i \omega \rho_0 q(\mathbf{r})$

where $k = \frac{2\pi f}{c}$ is the wavenumber, $c = 343\text{ m/s}$ is the speed of sound in air, and $\rho_0 = 1.204\text{ kg/m}^3$ is air density.

At room boundaries (rigid drywall, concrete, or glass), sound waves reflect almost completely. This means the normal particle velocity at the boundary face is zero. Mathematically, we impose homogeneous Neumann boundary conditions:

$\frac{\partial p(\mathbf{r})}{\partial n} = 0 \quad \text{on } \partial \Omega$

where $n$ is the unit normal vector pointing out of the room surface.

Step 1: Voxelizing Arbitrary 3D Geometries

Shoebox analytical equations ($f_{n_x, n_y, n_z} = \frac{c}{2}\sqrt{(n_x/L_x)^2 + (n_y/L_y)^2 + (n_z/L_z)^2}$) break down the moment a room has an L-extension, a dormer window, or an angled wall.

To support real-world room shapes, we discretize the 3D domain into a uniform voxel grid with a spatial step $h = 0.15\text{ m}$ by default (roughly $\lambda/11$ at the 200 Hz ceiling — the grid coarsens automatically for very large rooms to keep the eigenproblem tractable).

import numpy as np
from scipy.sparse import csr_matrix

def build_voxel_laplacian(air_voxels_mask, h=0.15):
    """
    Constructs the 7-point negative finite-difference Laplacian matrix
    with Neumann boundary conditions on arbitrary interior air voxels.
    """
    # Map (x, y, z) 3D grid index to a continuous 1D index
    grid_shape = air_voxels_mask.shape
    air_indices = np.argwhere(air_voxels_mask)
    voxel_to_idx = {tuple(coord): i for i, coord in enumerate(air_indices)}
    num_air = len(air_indices)
    
    row_ind, col_ind, data = [], [], []
    
    # 6 cardinal neighbor direction vectors
    neighbors = [(-1,0,0), (1,0,0), (0,-1,0), (0,1,0), (0,0,-1), (0,0,1)]
    
    for i, (x, y, z) in enumerate(air_indices):
        active_neighbors = 0
        for dx, dy, dz in neighbors:
            neighbor_coord = (x + dx, y + dy, z + dz)
            if neighbor_coord in voxel_to_idx:
                active_neighbors += 1
                j = voxel_to_idx[neighbor_coord]
                # Off-diagonal element for interior coupling
                row_ind.append(i)
                col_ind.append(j)
                data.append(-1.0 / (h**2))
        
        # Diagonal element enforces Neumann boundary implicitly
        # (ghost cell values equal boundary cell values)
        row_ind.append(i)
        col_ind.append(i)
        data.append(active_neighbors / (h**2))
        
    L = csr_matrix((data, (row_ind, col_ind)), shape=(num_air, num_air))
    return L, air_indices

By adjusting the diagonal term to match the exact number of active interior neighbors ($N_{\text{adj}} \le 6$), Neumann zero-flux boundary conditions are enforced naturally without constructing explicit ghost cells outside the mesh.

The resulting Laplacian matrix $\mathbf{L}$ is sparse, real, symmetric, and positive semi-definite.

Step 2: Extracting Mode Shapes with a Shift-and-Invert Krylov Solver

We don't need all the eigenvalues of the grid — the full grid can run into the hundreds of thousands of interior air cells. We only care about the modes below the 200 Hz handover to ray tracing.

The production solver runs on SLEPc's Krylov–Schur method with a MUMPS shift-invert factorization (via PETSc), not a plain SciPy call — MUMPS's sparse direct factorization is what makes shift-invert tractable at the grid sizes real rooms produce. The conceptual shape of the problem is the same one scipy.sparse.linalg.eigsh would solve on a small grid, so we show it here for illustration:

from scipy.sparse.linalg import eigsh

def eigensolve_room_modes(L, num_modes=60, c=343.0):
    """
    Illustrative small-grid version of the eigenproblem. Production runs
    SLEPc Krylov-Schur with MUMPS shift-invert (PETSc), which scales to
    the voxel counts real, unshoebox rooms produce.

    Solves L * phi = lambda * phi for the lowest spatial modes.
    Returns natural frequencies (Hz) and mode shape vectors.
    """
    # Shift-and-invert spectral transformation around sigma = 1e-5
    eigenvalues, eigenvectors = eigsh(
        L, 
        k=num_modes, 
        which='SM', 
        sigma=1e-5
    )
    
    # Convert spatial Laplacian eigenvalues (k^2) to frequency (Hz)
    # lambda = k^2 = (2 * pi * f / c)^2  ==>  f = (sqrt(lambda) * c) / (2 * pi)
    k_vals = np.sqrt(np.maximum(eigenvalues, 0))
    freqs_hz = (k_vals * c) / (2.0 * np.pi)
    
    return freqs_hz, eigenvectors

Each column of eigenvectors represents an unperturbed 3D mode shape $\phi_m(\mathbf{r})$. We normalize these vectors against the total room volume $V$ so that:

$\int_{\Omega} \phi_m^2(\mathbf{r}) \, dV = 1$

Step 3: First-Order Acoustic Perturbation for Absorbing Panels

What happens when a user drops a $50\text{ mm}$ thick porous absorber onto a side wall? Re-running the eigensolve for every single panel move would take seconds to minutes and destroy UI responsiveness.

Instead, we use First-Order Acoustic Modal Perturbation Theory, expressed as a per-mode quality factor $Q_n$:

$Q_n = \frac{4 \omega_n}{c \left( \sum \alpha \int |\phi_n|^2 \, dA + \text{panel terms} \right)}$

Panels contribute directly to the denominator as absorption-weighted surface integrals over the panel's placement, without touching the underlying eigenvectors. The updated modal decay time per third-octave band follows from $Q_n$:

$T_{60, n} = \frac{2.2 \cdot Q_n}{f_n}$

This update is cheap relative to the eigensolve itself, and the result is cached per room geometry — the same bare-room eigenpairs are reused across panel placements. The bare-room eigensolve itself is the expensive step: it can take anywhere from a couple of seconds to a couple of minutes depending on room size and shape, which is why it only runs when the geometry actually changes and is cached afterward. Once a geometry's modes are cached, sliding a panel around updates the low-frequency response in the UI without triggering a new solve.

Comparing Simulation Models across Regimes

Try the Interactive Simulation

We wrapped this entire solver pipeline into a fast, free interface at roomtreatment.diy. The 3D room render and standing-wave visualization run in WebGL in your browser; the physics itself — the voxel eigensolve and ray tracing — runs on our servers and streams results back to the viewer. You can input your room dimensions, drag your speakers and listening desk around, and inspect the 3D standing wave heatmaps and frequency response nulls for your exact space.

No installation required, no signup wall to run your first room calculation.

Hacker News Excerpt (Show HN)

Title: Show HN: roomtreatment.diy – Free web-based 3D room acoustics simulator

Post Body: Hey HN, I built roomtreatment.diy—a free web tool that simulates 3D room acoustics and standing waves.

Most web acoustic calculators use ray tracing or simple shoebox equations. But ray tracing completely breaks down in small rooms below 200 Hz, where sound acts as discrete standing waves rather than geometric rays.

To solve this for non-rectangular rooms (L-shapes, slanted ceilings), the backend discretizes the space into a 3D voxel grid ($h=0.15\text{ m}$ pitch) and solves the spatial negative Laplacian matrix $-\nabla^2$ under Neumann boundary conditions with a SLEPc Krylov-Schur eigensolver (MUMPS shift-invert). When you place wall absorbers or corner bass traps, it uses 1st-order modal perturbation theory to recompute modal $T_{60}$ decay rates and listener nulls without a full re-solve.

It's completely free, the 3D scene renders client-side in WebGL, and it requires no login to build a room model. I'd love your feedback on the physics model and UI experience!