# pffdtd: FDTD wave-based per l'acustica delle sale > Il fork ST-LINE del solver FDTD pffdtd di Brian Hamilton per l'acustica delle sale: stencil, CFL, bordi a impedenza, staircasing, costo computazionale e le ottimizzazioni CUDA a schema numerico invariato. Pubblicato: 2026-06-25 Aggiornato: 2026-08-25 Ambito: acustica Riferimento normativo: FDTD wave-based per l'acustica delle sale Repository: Pagina: --- 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](https://github.com/stefanofante/pffdtd); la scheda di vetrina è in [Open Lab](/open-lab/pffdtd/). ## 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, \qquad \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: | f_max | Δx | Nodi | Memoria | Passi | Aggiornamenti totali | |---|---|---|---|---|---| | 250 Hz | 0,2287 m | 4,2 · 10^4 | 1 MB | 3897 | 1,6 · 10^8 | | 500 Hz | 0,1143 m | 3,3 · 10^5 | 8 MB | 7794 | 2,6 · 10^9 | | 1000 Hz | 0,0572 m | 2,7 · 10^6 | 64 MB | 15588 | 4,2 · 10^10 | | 2000 Hz | 0,0286 m | 2,1 · 10^7 | 514 MB | 31176 | 6,7 · 10^11 | | 4000 Hz | 0,0143 m | 1,7 · 10^8 | 4,1 GB | 62353 | 1,1 · 10^13 | | 8000 Hz | 0,0071 m | 1,4 · 10^9 | 32,9 GB | 124707 | 1,7 · 10^14 | **Figura:** Grafico della memoria richiesta in funzione della frequenza massima Memoria della griglia in funzione della frequenza massima simulata, per una sala di 500 m³, in scala doppio-logaritmica. La retta ha pendenza 3: ogni ottava costa otto volte la memoria. La riga tratteggiata è il tetto di una macchina da 32 GB, raggiunto poco sotto gli 8 kHz. 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: | Volume | 8 GB | 32 GB | 128 GB | |---|---|---|---| | 50 m³ (regia) | 10,8 kHz | 17,1 kHz | 27,1 kHz | | 500 m³ (aula magna) | 5,0 kHz | 7,9 kHz | 12,6 kHz | | 15 000 m³ (sala da concerto) | 1,6 kHz | 2,6 kHz | 4,1 kHz | 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](/wiki/reverberation-time/)): 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](https://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](https://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](https://doi.org/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](https://doi.org/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](https://doi.org/10.1121/1.1915175). - **M. R. Schroeder**, sulla frequenza di transizione fra regime modale e diffuso. DOI [10.1121/1.1909343](https://doi.org/10.1121/1.1909343). - **L. Aspöck et al.**, *BRAS — Benchmark for Room Acoustical Simulation*, TU Berlin, 2020. DOI [10.14279/depositonce-6726.3](https://doi.org/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.*