<!-- canonical: https://jmrplens.github.io/phonometry/es/guides/fdtd-simulation/ -->
Source: https://jmrplens.github.io/phonometry/es/guides/fdtd-simulation/

La mayor parte de phonometry predice un **número**: un nivel, un tiempo de
reverberación, una pérdida por transmisión. El dominio `simulation` calcula el
**propio campo de ondas**: `fdtd_simulation` integra las ecuaciones acústicas
lineales sobre una malla 2D con el método de **diferencias finitas en el
dominio del tiempo (FDTD)**, de modo que la reflexión, la difracción, la
interferencia, la refracción en medios inhomogéneos y el comportamiento modal
emergen de primeros principios en lugar de modelarse término a término. La
implementación sigue la formulación de referencia para sonido en exteriores de
Attenborough y Van Renterghem, *Predicting Outdoor Sound* (2.ª ed., CRC Press
2021), capítulo 4: el esquema presión-velocidad escalonado en espacio y en
tiempo (Ecs. 4.11-4.12), la condición de estabilidad de Courant
(Ecs. 4.13-4.14), los contornos rígidos como velocidad normal nula en la cara
(Ec. 4.32) y el contorno de impedancia real independiente de la frecuencia
(Ecs. 4.33-4.35).

El solucionador es **determinista por diseño**: aritmética float64, sin
números aleatorios y avance numpy monohilo, de modo que las mismas entradas
producen salidas idénticas bit a bit en la misma plataforma. Es el motor de
las animaciones FDTD de esta documentación, promocionado a API pública con
fuentes, sondas de presión, obstáculos rasterizados, condiciones de contorno
por lado y un objeto de resultado congelado.

## 1. El esquema: una ecuación de onda sobre una malla

En un medio en reposo, las ecuaciones linealizadas de la dinámica de fluidos
se reducen a un sistema de primer orden en la presión acústica ``p`` y la
velocidad de partícula ``v`` (Attenborough y Van Renterghem, Ecs. 4.3-4.4):

$$
\frac{\partial p}{\partial t} + \rho c^2\,\nabla\!\cdot\!\mathbf{v} = 0,
\qquad
\frac{\partial \mathbf{v}}{\partial t} + \frac{1}{\rho}\,\nabla p = 0 .
$$

FDTD discretiza ambas sobre una **malla escalonada** (el análogo acústico de
la celda de Yee): la presión vive en los centros de celda y cada componente
de velocidad en las caras, a media celda, y los dos campos avanzan en
**leapfrog**, medio paso temporal desfasados (Ecs. 4.11-4.12). Evaluar cada
gradiente espacial exactamente donde lo necesita el otro campo cuadruplica la
precisión frente a una malla colocalizada (Ec. 4.9 frente a 4.10) y permite
actualizaciones in situ. Como solo se almacenan las caras interiores, el
borde del dominio es una **pared perfectamente rígida** (velocidad normal
nula, Ec. 4.32) salvo que se pida otro contorno.

El esquema explícito solo es estable mientras un frente de onda cruce como
mucho una celda por paso temporal. Con celdas cuadradas, el **número de
Courant** (Ec. 4.13) es

$$
\mathrm{CN} = c\,\Delta t\sqrt{\frac{1}{\Delta x^2} + \frac{1}{\Delta y^2}}
            = \frac{c\,\Delta t\,\sqrt{2}}{\Delta x} \le 1,
$$

y `fdtd_simulation` deriva el paso temporal del parámetro `cfl` (el número de
Courant, por defecto 0,6) y de la mayor velocidad del sonido del mapa; los
valores fuera de ``(0, 1)`` se rechazan porque el esquema carece de sentido
más allá del límite (Ec. 4.14).

```python
from phonometry import simulation

# Un dominio de aire de 3.0 x 2.0 m: 300 x 200 celdas de 1 cm.
res = simulation.fdtd_simulation(
    343.0, 0.01, 2.0e-3, shape=(200, 300),
    sources=[simulation.GaussianPulse(ix=60, iy=100, width=3.0e-4)],
    probes=[(200, 100)],
)
print(res.size)                  # (3.0, 2.0)  metros
print(round(res.dt * 1e6, 2))    # 12.37  microsegundos (CN = 0.6)
res.plot()                       # historias de presión en las sondas (figura del §3)
```

La malla es de índices: la celda ``(ix, iy)`` tiene su centro en
``((ix + 0.5) * dx, (iy + 0.5) * dx)`` metros, con las filas dibujadas hacia
abajo (la convención de ``imshow``), de modo que una posición en metros se
convierte con ``ix = round(x / dx - 0.5)``.

## 2. Fuentes, sondas, obstáculos y contornos

Tres tipos de fuente inyectan una **fuente blanda** (una contribución de
presión aditiva que no dispersa las ondas que pasan) en una celda de la
malla: `GaussianPulse` (un pulso de banda ancha con semianchura temporal
`width`), `CWSource` (un tono senoidal con rampa de coseno alzado para que su
arranque no salpique el campo con un transitorio de banda ancha) y
`SignalSource` (una forma de onda muestreada arbitraria, interpolada
linealmente sobre los pasos temporales de la simulación). Las sondas
registran la presión en su celda en cada paso temporal dentro del resultado.

La geometría se **rasteriza**: `obstacle_mask` marca celdas rígidas, y toda
cara en contacto con una celda marcada se cierra (de nuevo la Ec. 4.32), así
que paredes, barreras y difusores de cualquier forma son simples arrays
booleanos. Cada lado del dominio puede llevar su propia condición de
contorno:

- ``"rigid"`` (por defecto): un reflector perfecto, ``R = +1``.
- ``"absorbing"``: una capa esponja de `absorbing_layer_cells` celdas cuya
  tasa de absorción crece cuadráticamente, emulando un contorno abierto (el
  precursor sencillo de las capas perfectamente adaptadas de la sección
  4.2.3).
- una **impedancia específica real** ``Z`` en Pa·s/m (un escalar o un valor
  por celda de borde): el contorno de reacción local de las Ecs. 4.33-4.35,
  actualizado implícitamente, con coeficiente de reflexión a incidencia
  normal ``R = (Z - ρc)/(Z + ρc)``; ``Z = ρc`` es anecoico.

El motor de avance `FDTD2D` también es público: expone `step()`, `run()`,
los arrays del campo y la energía, para quien necesite acceso fotograma a
fotograma (las animaciones de la documentación lo usan directamente). Un
pulso plano lanzado por un conducto contra un borde de impedancia reproduce
el coeficiente de reflexión de libro:

```python

from phonometry import simulation

rho, c, dx = 1.2, 343.0, 0.01
sim = simulation.FDTD2D(c, dx, rho=rho, shape=(3, 1200),
                        edge_impedance={"right": 3.0 * rho * c})
x = (np.arange(1200) + 0.5) * dx
sim.p[:] = np.exp(-(((x - 6.0) / 0.15) ** 2))[None, :]   # pulso plano

trace = []
for _ in range(int(round(0.032 / sim.dt))):
    sim.step()
    trace.append(sim.p[1, 900])
trace = np.asarray(trace)
t = (np.arange(trace.size) + 1) * sim.dt
t_return = 6.0 / c + 3.0 / c        # vía la pared, de vuelta a x = 9 m
incident = trace[t < t_return - 0.001].max()
echo = trace[t > t_return]
print(round(float(echo[np.abs(echo).argmax()] / incident), 2))  # 0.5
# (Z - rho c)/(Z + rho c) = (3 - 1)/(3 + 1) = +0.5
```

## 3. Cuándo usarlo, y los límites del 2D

FDTD amortiza su coste cuando la **geometría gobierna la física**: difracción
en torno a una barrera o a través de una abertura, interferencia de caminos
directo y reflejado, dispersión por obstáculos, comportamiento modal de
recintos de forma irregular, refracción en un gradiente de velocidad del
sonido. Una sola ejecución captura **todas las frecuencias a la vez** (un
pulso excita toda la banda; una FFT de una sonda da el espectro), donde un
método en el dominio de la frecuencia necesita una resolución por frecuencia.

<details>
<summary>Mostrar el código de esta figura</summary>

```python

from phonometry import simulation

# Un campo libre de 3.0 x 2.0 m (bordes absorbentes) con una barrera
# rígida delgada: la sonda A ve el pulso directo más la reflexión de la
# barrera, la sonda B queda en sombra y solo recibe la onda difractada.
mask = np.zeros((200, 300), dtype=bool)
mask[60:, 150:154] = True
res = simulation.fdtd_simulation(
    343.0, 0.01, 9.0e-3, shape=(200, 300),
    sources=[simulation.GaussianPulse(ix=60, iy=100, width=3.0e-4)],
    probes=[(100, 100), (240, 100)],
    obstacle_mask=mask,
    boundaries="absorbing", absorbing_layer_cells=30,
    snapshot_every=75,
)

fig, (ax_f, ax_p) = plt.subplots(
    1, 2, figsize=(12.5, 5.0), gridspec_kw={"width_ratios": [1.25, 1.0]})
res.plot(kind="snapshot", frame=7, ax=ax_f, language="es")
res.plot(ax=ax_p, language="es")
plt.tight_layout()
plt.show()
```

</details>

A la inversa, cuando existe una forma cerrada validada (reverberación
estadística, atenuación en exteriores de ISO 9613-2, fuentes imagen en una
sala rectangular), la forma cerrada es miles de veces más barata: este
solucionador es el contraste y el demostrador, no el sustituto. El oráculo
funciona en ambos sentidos; una ejecución en caja rígida reproduce los modos
analíticos de la sala:

```python

from phonometry import simulation

lx, ly, dx = 1.0, 0.7, 0.02
nx, ny = round(lx / dx), round(ly / dx)
res = simulation.fdtd_simulation(
    343.0, dx, 0.35, shape=(ny, nx),
    sources=[simulation.GaussianPulse(ix=7, iy=5, width=2.0e-4)],
    probes=[(nx - 4, ny - 3)],
)
p = res.pressures[0]
spec = np.abs(np.fft.rfft(p * np.hanning(p.size), n=8 * p.size))
freqs = np.fft.rfftfreq(8 * p.size, res.dt)
sel = (freqs > 250) & (freqs < 350)
print(round(0.5 * 343.0 * float(np.hypot(1 / lx, 1 / ly)), 1))  # 299.1  modo (1,1) exacto
print(round(float(freqs[sel][np.argmax(spec[sel])]), 1))        # 298.9  medido
```

*El espectro en la sonda del ensayo de caja rígida presenta sus picos
exactamente sobre las frecuencias modales analíticas
$f = (c/2)\sqrt{(n_x/L_x)^2 + (n_y/L_y)^2}$ (Kuttruff 6.ª ed., cap. 3). El
desplazamiento apenas visible hacia la izquierda de los picos más altos es
la dispersión numérica del §4: las longitudes de onda cortas viajan
ligeramente lentas en la malla, así que las resonancias modeladas leen por
debajo en una fracción de porcentaje.*

<details>
<summary>Mostrar el código de esta figura</summary>

```python

from phonometry import simulation

lx, ly, dx, c = 1.0, 0.7, 0.02, 343.0
nx, ny = round(lx / dx), round(ly / dx)
res = simulation.fdtd_simulation(
    c, dx, 0.35, shape=(ny, nx),
    sources=[simulation.GaussianPulse(ix=7, iy=5, width=2.0e-4)],
    probes=[(nx - 4, ny - 3)],
)

# Una línea: la historia de presión en la sonda de la que sale el espectro.
res.plot(language="es")
plt.show()

# La comprobación modal: espectro en la sonda frente a los modos analíticos.
p = res.pressures[0]
spec = np.abs(np.fft.rfft(p * np.hanning(p.size), n=8 * p.size))
freqs = np.fft.rfftfreq(8 * p.size, res.dt)
sel = (freqs >= 100) & (freqs <= 450)
fig, ax = plt.subplots()
ax.plot(freqs[sel], 20 * np.log10(spec[sel] / spec[sel].max()))
for mx, my in [(1, 0), (0, 1), (1, 1), (2, 0), (2, 1)]:
    ax.axvline(0.5 * c * np.hypot(mx / lx, my / ly), ls=":", color="tab:red")
ax.set(xlabel="Frecuencia [Hz]", ylabel="Espectro en la sonda [dB re máx]",
       ylim=(-60, 6))
plt.show()
```

</details>

El dominio es **bidimensional**, y eso cambia la física, no solo el coste.
Una fuente puntual 2D es físicamente una **fuente lineal** infinita: su
amplitud se expande cilíndricamente como ``1/sqrt(r)`` (3,0 dB por
duplicación de distancia) en lugar del ``1/r`` esférico (6,0 dB) de una
fuente puntual 3D, y la respuesta impulsional 2D arrastra una estela tras el
frente en vez de pasar limpiamente. Los *patrones* de interferencia y
difracción son fieles; los niveles absolutos y las tasas de caída no son los
de una sala 3D. Trate las ejecuciones 2D como secciones y demostraciones, y
valide cualquier afirmación cuantitativa 3D contra una forma cerrada o un
solucionador 3D.

## 4. Dispersión numérica y precisión

La malla discreta propaga cada frecuencia a una velocidad ligeramente
errónea: las longitudes de onda cortas se retrasan, así que un pulso agudo
desarrolla una cola de rizado y las resonancias se desplazan un poco. Esta
**dispersión numérica** es la contrapartida discreta de la Ec. 4.15; sobre
los ejes de una malla cuadrada la relación de dispersión del esquema es

$$
\sin\!\left(\frac{\omega\,\Delta t}{2}\right)
  = \frac{c\,\Delta t}{\Delta x}\,
    \sin\!\left(\frac{k\,\Delta x}{2}\right),
$$

con un error relativo de frecuencia a primer orden de magnitud
``(1 - S^2) (k dx)^2 / 24`` sobre los ejes de la malla (la frecuencia
modelada queda por debajo de la real, así que el error con signo es
negativo), con ``S = c dt / dx``; el error es máximo exactamente sobre un
eje y se anula a lo largo de la diagonal de las celdas en el límite de
Courant ``CN = 1``. La regla práctica es resolver **al menos 10 celdas por
longitud de onda más corta**, ``dx <= c_min / (10 f_max)`` con la velocidad
de sonido más baja del dominio: con exactamente 10 celdas la cota de
Courant pequeño ``(k dx)^2 / 24`` da alrededor del 1,6 %, que el factor
``1 - S^2`` reduce a en torno al 1,4 % con el ``cfl = 0.6`` por defecto
(en un dominio heterogéneo el paso temporal lo fijan las celdas más
rápidas, así que las regiones lentas trabajan a un número de Courant local
menor y quedan más cerca de la cota del 1,6 %), y toda componente mejor
resuelta o fuera de eje es más precisa. Con
``dx = 1`` cm el punto de 10 celdas en aire queda hacia los 3,4 kHz, y
reducir ``dx`` a la mitad divide el error entre cuatro (el esquema es de
segundo orden, y la batería de validación mide ese orden observado bajo
refinamiento de malla). Los tests anclan el solucionador a oráculos
analíticos: autofrecuencias de caja y de conducto, tiempos de llegada en
campo libre y caída cilíndrica, el eco de la fuente imagen de una pared
rígida, el coeficiente de reflexión de impedancia anterior y la propia
relación de dispersión.

El `FDTDResult` congelado lleva el eje temporal, las historias de presión por
sonda, las posiciones de las sondas en metros, los metadatos de la malla, las
fuentes, las instantáneas opcionales del campo con sus tiempos y la máscara
de obstáculos; su `.plot()` dibuja las historias de las sondas, y
`.plot(kind="snapshot")` representa un campo registrado con la geometría
superpuesta.

## Qué cubre esta guía

**Cubierto.** El esquema FDTD presión-velocidad escalonado de Attenborough y
Van Renterghem, capítulo 4: las ecuaciones de gobierno (Ecs. 4.3-4.4), la
actualización leapfrog (Ecs. 4.11-4.12), la condición de estabilidad de
Courant (Ecs. 4.13-4.14), el contorno rígido (Ec. 4.32) y el contorno de
impedancia real independiente de la frecuencia (Ecs. 4.33-4.35), expuestos a
través de `fdtd_simulation`/`FDTD2D` con las fuentes `GaussianPulse`,
`CWSource` y `SignalSource`, sondas de presión y una `obstacle_mask`
rasterizada. Validado contra oráculos en forma cerrada: los modos normales de
sala rígida (Kuttruff cap. 3), los tiempos de llegada en campo libre y el
decaimiento cilíndrico en 2D, el eco de imagen de pared rígida, el
coeficiente de reflexión de impedancia y la relación de dispersión del §4,
con el orden de convergencia medido coincidiendo con el diseño de segundo
orden del esquema.

**No cubierto.** El solver es solo bidimensional: modela una fuente lineal en
sección transversal (propagación cilíndrica `1/sqrt(r)`), no la propagación
esférica `1/r` de una fuente puntual 3D, así que los niveles absolutos y las
tasas de decaimiento no son los de una sala 3D. El contorno abierto es la
capa absorbente de rampa cuadrática descrita como «el precursor sencillo» de
una verdadera capa perfectamente adaptada, no una PML en sí. Las ecuaciones
de gobierno asumen un medio en reposo, así que no se modela la advección por
viento o flujo, y el único contorno de impedancia es el real independiente
de la frecuencia de las Ecs. 4.33-4.35.

## Véase también

- [Informe de conformidad](/phonometry/es/reference/conformance/): dos de los anclajes
  de validación en forma cerrada del §4 se ejecutan ahí.
- Referencia de la API: [`simulation.fdtd`](/phonometry/es/reference/api/simulation/fdtd/).

## Respuestas rápidas

### ¿Cómo elijo el paso de malla en FDTD?

Resuelva al menos 10 celdas por longitud de onda más corta,
``dx <= c_min / (10 f_max)``, con la velocidad de sonido más baja del
dominio. Con exactamente 10 celdas la cota del error de dispersión sobre el
eje es alrededor del 1,6 %, reducida a en torno al 1,4 % con el
``cfl = 0.6`` por defecto, y reducir ``dx`` a la mitad divide el error entre
cuatro (el esquema es de segundo orden). Con ``dx = 1`` cm el punto de 10
celdas en aire queda hacia los 3,4 kHz.

### ¿Qué número de Courant mantiene estable una simulación FDTD?

El esquema explícito solo es estable mientras un frente de onda cruce como
mucho una celda por paso temporal. Con celdas cuadradas el número de Courant
es $\mathrm{CN} = c\,\Delta t\,\sqrt{2}/\Delta x \le 1$ (Attenborough y Van
Renterghem, Ec. 4.13). `fdtd_simulation` deriva el paso temporal del
parámetro `cfl` (el número de Courant, por defecto 0,6) y de la mayor
velocidad del sonido del mapa, y rechaza los valores fuera de $(0, 1)$.

### ¿Puedo fiarme de los niveles absolutos de una ejecución FDTD 2D?

No. Una fuente puntual 2D es físicamente una fuente lineal infinita: su
amplitud se expande cilíndricamente como ``1/sqrt(r)``, 3,0 dB por
duplicación de distancia, en lugar del ``1/r`` esférico (6,0 dB por
duplicación) de una fuente puntual 3D. Los patrones de interferencia y
difracción son fieles, pero los niveles absolutos y las tasas de caída no
son los de una sala 3D; valide cualquier afirmación cuantitativa 3D contra
una forma cerrada o un solucionador 3D.