pffdtd: FDTD wave-based per l'acustica delle sale

Questa è la nota tecnica del fork pffdtd che ST-LINE mantiene: cosa fa un solver FDTD wave-based per l’acustica delle sale, quanto costa, dove sbaglia, e che cosa il fork ha cambiato rispetto all’originale di Brian Hamilton — lo schema numerico no, il motore CUDA sì. Il repository è su github.com/stefanofante/pffdtd; la scheda di vetrina è in Open Lab.

L’FDTD per l’acustica delle sale

Il Finite-Difference Time-Domain (FDTD) risolve l’acustica lineare nel dominio del tempo discretizzando direttamente le equazioni del primo ordine che legano pressione e velocità delle particelle. In assenza di flusso medio e per piccole perturbazioni, la conservazione della quantità di moto e della massa danno il sistema

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

dove $p$ è la pressione acustica, $\mathbf{v}$ la velocità delle particelle, $\rho$ la densità dell’aria e $c$ la velocità del suono. Eliminando la velocità si ottiene l’equazione delle onde scalare ${\partial }_{tt}p=c^2{\nabla }^{2}p$, ed è questa che l’FDTD avanza passo dopo passo su una griglia regolare, sostituendo le derivate con differenze finite centrate nel tempo e nello spazio.

Stencil cartesiano vs FCC

La scelta dello stencil — l’insieme di nodi vicini che approssimano il laplaciano — governa l’accuratezza e il costo. Lo schema più semplice è lo stencil a 7 punti su griglia cartesiana (il nodo centrale più i sei primi vicini lungo gli assi). È economico ma fortemente anisotropo: la velocità numerica delle onde dipende dalla direzione di propagazione rispetto alla griglia, e questa dispersione numerica introduce errore di fase che cresce con la frequenza.

Una griglia face-centred cubic (FCC) con uno stencil a 13 punti distribuisce i vicini in modo più isotropo. A parità di dispersione tollerata, l’FCC permette un passo spaziale più grosso e arriva a costare circa cinque volte meno memoria dello stencil cartesiano a 7 punti — un fattore decisivo, perché in FDTD il numero di nodi cresce con il cubo della frequenza massima e il consumo di RAM è il vincolo dominante.

Stabilità, CFL e conservazione dell’energia

Uno schema esplicito nel tempo è stabile solo se il passo temporale rispetta la condizione di Courant–Friedrichs–Lewy (CFL): l’informazione non può attraversare più di una cella per passo. In forma compatta,

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

con ${\lambda }_{\max }$ che dipende dallo stencil (è $1/\sqrt{3}$ per il 7 punti cartesiano in 3D). Superare il limite fa esplodere la soluzione; restare appena sotto massimizza l’accuratezza per nodo.

Una proprietà più profonda è la conservazione dell’energia. Lo schema può essere costruito in forma passiva, cioè con un funzionale di energia discreta che non cresce mai in assenza di sorgenti e dissipazione. In doppia precisione questa conservazione è verificabile a precisione di macchina: il bilancio energetico chiude a meno dell’epsilon numerico passo dopo passo, e questo è il controllo di sanità più severo che un solver wave-based possa offrire. In singola precisione la conservazione vale solo entro un margine più ampio, per cui l’engine adotta dei safeguard (monitoraggio dell’energia, scelte di ordine delle operazioni) per evitare derive lente nei run lunghi.

Bordi a impedenza frequenza-dipendente (ADE)

I materiali reali non hanno impedenza costante: assorbono in modo diverso a frequenze diverse. Modellare un bordo come condizione di impedenza $Z(\omega )$ dipendente dalla frequenza significa, nel dominio del tempo, una convoluzione fra pressione e risposta del bordo — costosa e con memoria lunga. La tecnica delle auxiliary differential equation (ADE) aggira la convoluzione: rappresenta l’impedenza come una funzione razionale di $\omega$ e introduce, per ogni nodo di bordo, alcune variabili ausiliarie che evolvono secondo una ODE locale. La convoluzione globale diventa così un piccolo sistema di ODE aggiornate in loco, compatibile con la natura esplicita e locale dell’FDTD.

Quanto costa l’FDTD, in numeri

Che il numero di nodi cresca col cubo della frequenza e che la RAM sia il vincolo dominante è vero, ma finché resta qualitativo non dice quali problemi siano affrontabili. Mettiamolo in cifre su una sala di 500 m³, con sei punti per lunghezza d’onda alla frequenza massima, 1,5 s di coda simulata, stencil cartesiano a 7 punti (CFL 1/√3) e tre campi in doppia precisione, cioè 24 byte per nodo:

La legge è semplice e brutale. Il passo spaziale scala come 1/f, quindi i nodi come f³; ma il CFL lega il passo temporale a quello spaziale, quindi anche i passi temporali crescono come f. Il lavoro totale va dunque come f⁴: salire di un’ottava costa otto volte la memoria e sedici volte il tempo. Fra 250 Hz e 8 kHz, cinque ottave, ci sono sei ordini di grandezza di memoria e sei di lavoro.

È questo che rende l’FDTD un metodo per la parte bassa dello spettro, non una scelta stilistica.

Cosa compra lo stencil FCC

Detto così, il fattore cinque di memoria dell’FCC sembra un’ottimizzazione fra le tante. In realtà, dentro una legge cubica, cinque volte la memoria vale 5^(1/3) = 1,71 volte la frequenza massima, cioè circa tre quarti di ottava di banda in più a parità di macchina. E poiché un passo spaziale 1,71 volte più grosso porta con sé anche 1,71 volte meno passi temporali, il risparmio sul lavoro totale è di 8,5 volte, non di cinque.

Tre quarti di ottava sono la differenza fra fermarsi a 5 kHz e arrivare a 8,5 kHz, oppure fra un run che sta in 32 GB e uno che ne chiede 160. Detto in questi termini si capisce perché la scelta dello stencil venga prima di ogni micro-ottimizzazione del kernel.

Dove sta il soffitto, in pratica

Il tetto di banda dipende dal volume, e la dipendenza è cubica anche lì. Con le stesse ipotesi, la frequenza massima che entra in memoria:

Va confrontato con la frequenza sotto la quale un metodo a onde è davvero necessario, cioè la frequenza di Schroeder, dove i modi propri si separano e la descrizione statistica cade (vedi la wiki sul tempo di riverberazione): 179 Hz per la regia da 50 m³ con T60 0,4 s, 89 Hz per l’aula con T60 1 s, 23 Hz per la sala da concerto con T60 2 s.

Il confronto è istruttivo, e va nel verso opposto a quello che si aspetterebbe: il soffitto dell’FDTD sta una o due decadi sopra la frequenza in cui un metodo a onde diventa indispensabile. Su una macchina ordinaria si copre a onde molto più spettro di quanto la fisica del campo diffuso richieda.

Da qui la vera ragione degli stack ibridi, che non è l’impossibilità di arrivare a f_S. È che (a) si vuole accuratezza a onde ben sopra f_S, dove contano le prime riflessioni e la diffrazione, e la geometria acustica è approssimativa proprio lì; e (b) nel design inverso la simulazione forward non si esegue una volta ma centinaia. Un run a 2 kHz sulla sala da 500 m³ costa 6,7 · 10¹¹ aggiornamenti di nodo: un’ottimizzazione con cento iterazioni e due solve per gradiente ne chiede 1,3 · 10¹⁴. È in quel regime che il costo per grado di libertà diventa il parametro che decide tutto — ed è il motivo per cui la sezione seguente arriva al Galerkin discontinuo.

Lo staircasing e la sua correzione

Una griglia regolare rappresenta bene pareti allineate agli assi, ma una superficie inclinata o curva viene approssimata “a scalini” (staircasing). Il problema non è solo estetico: la superficie di bordo effettiva vista dallo schema è sistematicamente sbagliata — gli scalini gonfiano l’area esposta — e poiché l’assorbimento totale è proporzionale all’area, lo staircasing porta a sovrastimare l’assorbimento e quindi, a volte, a stimare male i tempi di decadimento. Su geometrie sfavorevoli l’errore sul tempo di riverbero può arrivare intorno al 50 %.

La correzione per area effettiva rimedia pesando il contributo di ogni nodo di bordo con l’area realmente intercettata dalla superficie continua, anziché con l’area dello scalino. Applicata su griglie via via più fini, riporta l’errore sotto l’1 %. Resta però un rimedio: la geometria sottostante è ancora staircased, e questo è uno dei motivi strutturali per cui, quando serve la geometria esatta, si passa a un metodo body-conforming (vedi più sotto il Galerkin discontinuo).

Il fork pffdtd: cosa è stato ottimizzato

pffdtd nasce come simulatore FDTD di Brian Hamilton (University of Edinburgh, 2021), distribuito con licenza MIT. Il nostro fork (github.com/stefanofante/pffdtd) lascia lo schema numerico invariato: a parità di input, l’output forward del fork coincide con l’engine Python di riferimento a precisione di macchina. Non abbiamo toccato la fisica; abbiamo lavorato sotto, sul motore GPU CUDA, dove l’implementazione originale lasciava prestazioni e portabilità sul tavolo.

Le ottimizzazioni, in sintesi:

• Build per architettura nativa. Invece di compilare per sm_35 (Kepler) e affidarsi alla ricompilazione PTX-JIT a runtime sulle GPU moderne, il fork compila per l’architettura reale della scheda. Si eliminano la latenza di JIT al primo lancio e le inefficienze del codice generato per un’ISA obsoleta.
• Peer-access multi-GPU reale. L’originale, su problemi che eccedono una scheda, faceva staging silenzioso dei dati attraverso la RAM dell’host, con un collo di bottiglia sul bus PCIe. Il fork abilita il peer-access diretto fra GPU, quando il topology lo consente, così le partizioni del dominio comunicano device-to-device.
• Iniezione delle sorgenti in batch. L’originale lanciava un kernel <<<1,1>>> — un singolo thread — per ogni sorgente e per ogni timestep, sprecando il parallelismo della GPU. Il fork inietta le sorgenti in batch, in un kernel solo.
• Read-out a blocchi. L’estrazione dei ricevitori avviene a blocchi anziché punto-per-punto, riducendo il numero di trasferimenti.
• Budget di memoria interrogato a runtime. Invece di una soglia hardcodata, l’engine chiede alla GPU quanta memoria è disponibile (cudaMemGetInfo) e dimensiona il problema di conseguenza — scopre l’hardware e si adatta, invece di assumere.

Due profili hardware di tuning

Il tuning è stato fatto su due macchine deliberatamente diverse, una RTX 4500 Ada con 24 GB di VRAM dedicata e una DGX Spark GB10 con 128 GB di memoria unificata. La differenza nel comportamento in memoria è istruttiva e ha guidato la logica di dimensionamento. Su VRAM dedicata l’over-allocation fallisce duro: superata la capacità, l’allocazione va in errore e il run si interrompe. Su memoria unificata il sistema degrada gradualmente, paginando fra GPU e host: il run continua, ma rallenta. Interrogare il budget a runtime e dimensionare il dominio sotto la soglia reale evita entrambe le patologie — il crash secco da una parte, il degrado silenzioso dall’altra.

Riferimenti

• B. Hamilton, pffdtd — solver FDTD per acustica delle sale, 2021, licenza MIT. github.com/bsxfun/pffdtd.
• F. Mondet et al., modello di impedenza frazionaria a pochi parametri per superfici acustiche, 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, sulla frequenza di transizione fra regime modale e diffuso. 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.

I valori riportati provengono dalla documentazione di progetto e dal README del fork; i DOI non verificati direttamente sono citati per autore, anno e titolo senza DOI.