Agent Skills
› synthetic-sciences/openscience
› wave-propagation
wave-propagation
GitHub模拟声、电磁、弹性及量子波传播。支持1D/2D/3D波动方程，采用FDTD和谱方法，涵盖源类型、散射、色散分析及多种边界条件处理。
Trigger Scenarios
模拟声波、光波、地震波或水波传播
分析障碍物引起的波散射现象
进行共振和驻波分析
求解波导和腔体问题
时域脉冲传播模拟
Install
npx skills add synthetic-sciences/openscience --skill wave-propagation -g -y
SKILL.md
Frontmatter
{
"name": "wave-propagation",
"tags": [
"Wave Equation",
"FDTD",
"Acoustics",
"Electromagnetics",
"Propagation",
"Simulation"
],
"author": "Synthetic Sciences",
"license": "MIT",
"version": "1.0.0",
"category": "physics",
"description": "Simulate wave propagation — acoustic, electromagnetic, elastic, and quantum waves. FDTD, spectral methods, and absorbing boundary conditions for 1D\/2D\/3D wave equations with sources, scattering, and dispersion.",
"dependencies": [
"scipy>=1.11.0",
"numpy>=1.24.0",
"matplotlib>=3.7.0"
]
}
Wave Propagation
Overview
Simulate wave propagation using finite-difference time-domain (FDTD) and spectral methods. Supports acoustic, electromagnetic, and elastic waves with various source types, boundary conditions, and media.
When to Use
- Simulating sound, light, seismic, or water waves
- Wave scattering from obstacles
- Resonance and standing wave analysis
- Waveguide and cavity problems
- Time-domain pulse propagation
Core Workflows
1. 1D Wave Equation (FDTD)
$$\frac{\partial^2 u}{\partial t^2} = c^2 \frac{\partial^2 u}{\partial x^2}$$
import numpy as np
import matplotlib.pyplot as plt
def wave_1d_fdtd(Nx=1000, Nt=2000, c=1.0, dx=0.01, CFL=0.9,
source_pos=0.2, source_freq=5.0, bc='absorbing'):
"""
1D wave equation solver using FDTD (leapfrog in time).
CFL condition: c*dt/dx ≤ 1
"""
dt = CFL * dx / c
x = np.arange(Nx) * dx
# Fields: u at three time levels
u_prev = np.zeros(Nx)
u_curr = np.zeros(Nx)
u_next = np.zeros(Nx)
r2 = (c * dt / dx)**2 # CFL number squared
print(f"CFL = {c*dt/dx:.4f}, dt = {dt:.6f}")
source_idx = int(source_pos / dx)
snapshots = []
for n in range(Nt):
t = n * dt
# Interior update
u_next[1:-1] = (2*u_curr[1:-1] - u_prev[1:-1] +
r2 * (u_curr[2:] - 2*u_curr[1:-1] + u_curr[:-2]))
# Source: Ricker wavelet
t0 = 1.0 / source_freq
source = (1 - 2*(np.pi*source_freq*(t - t0))**2) * \
np.exp(-(np.pi*source_freq*(t - t0))**2)
u_next[source_idx] += dt**2 * source
# Boundary conditions
if bc == 'absorbing':
# Mur first-order ABC
u_next[0] = u_curr[1] + (c*dt - dx)/(c*dt + dx) * (u_next[1] - u_curr[0])
u_next[-1] = u_curr[-2] + (c*dt - dx)/(c*dt + dx) * (u_next[-2] - u_curr[-1])
elif bc == 'fixed':
u_next[0] = 0
u_next[-1] = 0
elif bc == 'periodic':
u_next[0] = u_next[-2]
u_next[-1] = u_next[1]
u_prev = u_curr.copy()
u_curr = u_next.copy()
if n % (Nt // 10) == 0:
snapshots.append((t, u_curr.copy()))
return x, snapshots
x, snaps = wave_1d_fdtd()
fig, ax = plt.subplots(figsize=(12, 6))
for t, u in snaps:
ax.plot(x, u + t*0.3, linewidth=0.8, label=f't={t:.3f}')
ax.set_xlabel('x [m]')
ax.set_ylabel('u (offset by time)')
ax.set_title('1D Wave Propagation (Ricker Source, Absorbing BCs)')
ax.legend(ncol=2, fontsize=8)
ax.grid(True, alpha=0.3)
plt.savefig('wave_1d.png', dpi=150, bbox_inches='tight')
2. 2D Wave Equation
def wave_2d_fdtd(Nx=200, Ny=200, Nt=500, c=1.0, dx=0.01, CFL=0.7):
"""2D wave equation using FDTD."""
dt = CFL * dx / (c * np.sqrt(2)) # 2D CFL
u_prev = np.zeros((Nx, Ny))
u_curr = np.zeros((Nx, Ny))
u_next = np.zeros((Nx, Ny))
r2 = (c * dt / dx)**2
source_i, source_j = Nx // 4, Ny // 2
snapshots = []
for n in range(Nt):
t = n * dt
# Interior
u_next[1:-1, 1:-1] = (
2*u_curr[1:-1, 1:-1] - u_prev[1:-1, 1:-1] + r2 * (
u_curr[2:, 1:-1] + u_curr[:-2, 1:-1] +
u_curr[1:-1, 2:] + u_curr[1:-1, :-2] -
4*u_curr[1:-1, 1:-1]
)
)
# Source
freq = 10.0
t0 = 0.1
source = (1 - 2*(np.pi*freq*(t-t0))**2) * np.exp(-(np.pi*freq*(t-t0))**2)
u_next[source_i, source_j] += dt**2 * source
# Absorbing BCs (simple)
u_next[0, :] = u_next[1, :]
u_next[-1, :] = u_next[-2, :]
u_next[:, 0] = u_next[:, 1]
u_next[:, -1] = u_next[:, -2]
u_prev = u_curr.copy()
u_curr = u_next.copy()
if n % (Nt // 6) == 0:
snapshots.append((t, u_curr.copy()))
return snapshots
snaps_2d = wave_2d_fdtd()
fig, axes = plt.subplots(2, 3, figsize=(15, 10))
for ax, (t, u) in zip(axes.flat, snaps_2d):
im = ax.imshow(u.T, cmap='RdBu_r', vmin=-0.01, vmax=0.01,
origin='lower', extent=[0, 2, 0, 2])
ax.set_title(f't = {t:.3f} s')
ax.set_xlabel('x [m]')
ax.set_ylabel('y [m]')
plt.suptitle('2D Wave Propagation')
plt.tight_layout()
plt.savefig('wave_2d.png', dpi=150, bbox_inches='tight')
3. Spectral Method (Periodic Waves)
def wave_spectral(N=256, L=2*np.pi, T=10, c=1.0, dt=0.01):
"""Wave equation via pseudospectral method (periodic BC)."""
x = np.linspace(0, L, N, endpoint=False)
k = np.fft.fftfreq(N, d=L/N) * 2 * np.pi
# Initial: Gaussian pulse
u = np.exp(-((x - L/2)**2) / 0.05)
v = np.zeros(N) # du/dt = 0
u_hat = np.fft.fft(u)
v_hat = np.fft.fft(v)
omega2 = (c * k)**2
Nt = int(T / dt)
snaps = [(0, u.copy())]
for n in range(Nt):
# Leapfrog in Fourier space
v_hat -= dt * omega2 * u_hat
u_hat += dt * v_hat
if (n+1) % (Nt // 8) == 0:
snaps.append(((n+1)*dt, np.real(np.fft.ifft(u_hat))))
return x, snaps
4. Dispersion Relation Analysis
def measure_dispersion(simulation_data, dx, dt):
"""
Measure dispersion relation from simulation data.
Compute 2D FFT (space-time) to get ω(k).
"""
# 2D FFT
ft = np.fft.fft2(simulation_data)
ft_shifted = np.fft.fftshift(ft)
Nt, Nx = simulation_data.shape
k = np.fft.fftshift(np.fft.fftfreq(Nx, d=dx)) * 2 * np.pi
omega = np.fft.fftshift(np.fft.fftfreq(Nt, d=dt)) * 2 * np.pi
plt.figure(figsize=(8, 6))
plt.pcolormesh(k, omega, np.log10(np.abs(ft_shifted)**2 + 1e-20),
cmap='hot', shading='auto')
plt.plot(k, np.abs(k), 'w--', linewidth=1, label='ω = c|k| (exact)')
plt.xlabel('Wavenumber k [rad/m]')
plt.ylabel('Frequency ω [rad/s]')
plt.title('Dispersion Relation')
plt.colorbar(label='log₁₀|FFT|²')
plt.legend()
plt.savefig('dispersion.png', dpi=150)
CFL Stability Conditions
| Dimension | Condition | Notes |
|---|---|---|
| 1D | c·dt/dx ≤ 1 | Exact for FDTD |
| 2D | c·dt/dx ≤ 1/√2 | Square grid |
| 3D | c·dt/dx ≤ 1/√3 | Cubic grid |
Source Types
| Source | Formula | Use For |
|---|---|---|
| Ricker wavelet | (1-2(πft₀)²)exp(-(πft₀)²) | Seismic, broadband pulse |
| Gaussian pulse | exp(-t²/2σ²) | Simple test pulse |
| Sine burst | sin(2πft) × window | Narrowband excitation |
| Point source | δ(x-x₀)·s(t) | Monopole radiation |
Troubleshooting
| Symptom | Fix |
|---|---|
| Solution blows up | CFL violation — reduce dt |
| Reflections from boundary | Use absorbing BC (Mur, PML) |
| Numerical dispersion | Reduce dx (need ~10-20 points per wavelength) |
| Spectral ringing | Smooth the source function |
Version History
- e9844a4 Current 2026-07-11 17:33
Dependencies
-
required
scipy>=1.11.0 -
required
numpy>=1.24.0 -
required
matplotlib>=3.7.0
