|
| 1 | +## 2D FDTD Implementation {#sec-em-2d-devito} |
| 2 | + |
| 3 | +Extending to two dimensions introduces new challenges: the 2D Yee cell has |
| 4 | +three field components arranged on cell edges and centers, and the stability |
| 5 | +condition becomes more restrictive. Most importantly, we need absorbing |
| 6 | +boundary conditions that work at oblique incidence angles. |
| 7 | + |
| 8 | +### TE and TM Modes in 2D |
| 9 | + |
| 10 | +In 2D simulations, we have two independent polarizations: |
| 11 | + |
| 12 | +**TE Mode** (Transverse Electric): $E_x$, $E_y$, $H_z$ |
| 13 | + |
| 14 | +- Electric field lies in the $xy$-plane |
| 15 | +- Magnetic field is normal to the plane |
| 16 | + |
| 17 | +**TM Mode** (Transverse Magnetic): $H_x$, $H_y$, $E_z$ |
| 18 | + |
| 19 | +- Magnetic field lies in the $xy$-plane |
| 20 | +- Electric field is normal to the plane |
| 21 | + |
| 22 | +We implement the TM mode (a single out-of-plane electric component). The governing equations are: |
| 23 | +$$ |
| 24 | +\frac{\partial H_x}{\partial t} = -\frac{1}{\mu}\frac{\partial E_z}{\partial y} |
| 25 | +$$ {#eq-em-2d-Hx} |
| 26 | +$$ |
| 27 | +\frac{\partial H_y}{\partial t} = \frac{1}{\mu}\frac{\partial E_z}{\partial x} |
| 28 | +$$ {#eq-em-2d-Hy} |
| 29 | +$$ |
| 30 | +\frac{\partial E_z}{\partial t} = \frac{1}{\varepsilon}\left(\frac{\partial H_y}{\partial x} - \frac{\partial H_x}{\partial y}\right) |
| 31 | +$$ {#eq-em-2d-Ez} |
| 32 | +
|
| 33 | +### The 2D Solver |
| 34 | +
|
| 35 | +The `src.em.maxwell2D_devito` module provides a 2D TM-mode solver |
| 36 | +using Devito [@devito-api]: |
| 37 | +
|
| 38 | +```python |
| 39 | +from src.em import solve_maxwell_2d, gaussian_source_2d |
| 40 | +import numpy as np |
| 41 | +
|
| 42 | +# Gaussian initial condition at domain center |
| 43 | +def E_init(X, Y): |
| 44 | + return gaussian_source_2d(X, Y, x0=0.5, y0=0.5, sigma=0.05) |
| 45 | +
|
| 46 | +result = solve_maxwell_2d( |
| 47 | + Lx=1.0, Ly=1.0, # Domain size [m] |
| 48 | + Nx=100, Ny=100, # Grid points |
| 49 | + T=3e-9, # Simulation time [s] |
| 50 | + CFL=0.5, # Courant number (<= 1/sqrt(2) for stability) |
| 51 | + E_init=E_init, |
| 52 | + pml_width=15, # PML absorbing boundary (grid cells) |
| 53 | + save_history=True, |
| 54 | +) |
| 55 | +``` |
| 56 | +
|
| 57 | +### 2D CFL Condition |
| 58 | +
|
| 59 | +The stability limit in 2D is more restrictive than 1D: |
| 60 | +$$ |
| 61 | +c \Delta t \leq \frac{1}{\sqrt{\frac{1}{\Delta x^2} + \frac{1}{\Delta y^2}}} |
| 62 | +$$ {#eq-em-cfl-2d-detailed} |
| 63 | +
|
| 64 | +For a uniform grid ($\Delta x = \Delta y = h$): |
| 65 | +$$ |
| 66 | +c \Delta t \leq \frac{h}{\sqrt{2}} \approx 0.707 h |
| 67 | +$$ |
| 68 | +
|
| 69 | +This means the maximum Courant number is $1/\sqrt{2} \approx 0.707$, compared |
| 70 | +to 1.0 in 1D. The solver enforces this: |
| 71 | +
|
| 72 | +```python |
| 73 | +>>> from src.em import solve_maxwell_2d |
| 74 | +>>> solve_maxwell_2d(Lx=1.0, Ly=1.0, Nx=100, Ny=100, T=1e-9, CFL=0.9) |
| 75 | +ValueError: CFL=0.9 > 1/sqrt(2) ~ 0.707 violates 2D stability condition |
| 76 | +``` |
| 77 | +
|
| 78 | +### Perfectly Matched Layer (PML) {#sec-em-pml} |
| 79 | +
|
| 80 | +Simple ABCs like Mur's condition work poorly in 2D because waves approach |
| 81 | +boundaries at various angles. The **Perfectly Matched Layer** (PML), |
| 82 | +introduced by Berenger in 1994 [@berenger1994], provides a much more |
| 83 | +effective solution. |
| 84 | +
|
| 85 | +The key insight is to create an absorbing region with: |
| 86 | +
|
| 87 | +1. **Matched impedance**: No reflection at the interface (for any angle) |
| 88 | +2. **Exponential decay**: Waves attenuate as they propagate into the PML |
| 89 | +
|
| 90 | +Mathematically, this is achieved through **complex coordinate stretching**: |
| 91 | +$$ |
| 92 | +\tilde{x} = x + \frac{j}{\omega}\int_0^x \sigma_x(x') dx' |
| 93 | +$$ |
| 94 | +
|
| 95 | +This makes the wave "see" a lossy medium while maintaining impedance matching. |
| 96 | +For the acoustic wave equation, an efficient PML formulation is given |
| 97 | +by @grote_sim2010. The Convolutional PML (CPML) variant [@roden_gedney2000] |
| 98 | +offers improved stability for anisotropic and dispersive media. |
| 99 | +
|
| 100 | +### PML Implementation |
| 101 | +
|
| 102 | +Our solver supports two PML implementations, selectable via the |
| 103 | +`pml_type` parameter: |
| 104 | +
|
| 105 | +**Graded-conductivity absorbing layer** (`pml_type='conductivity'`). |
| 106 | +The simplest approach adds a spatially varying conductivity to the |
| 107 | +update equations. The conductivity increases polynomially from zero |
| 108 | +at the PML interface to $\sigma_{\max}$ at the outer boundary: |
| 109 | +
|
| 110 | +```python |
| 111 | +def create_pml_profile(N, pml_width, sigma_max, order=3): |
| 112 | + """Create polynomial PML conductivity profile.""" |
| 113 | + sigma = np.zeros(N) |
| 114 | + for i in range(pml_width): |
| 115 | + d = (pml_width - i) / pml_width |
| 116 | + sigma[i] = sigma_max * (d ** order) # Left PML |
| 117 | + for i in range(N - pml_width, N): |
| 118 | + d = (i - (N - pml_width - 1)) / pml_width |
| 119 | + sigma[i] = sigma_max * (d ** order) # Right PML |
| 120 | + return sigma |
| 121 | +``` |
| 122 | +
|
| 123 | +**Convolutional PML** (`pml_type='cpml'`, default). The CPML |
| 124 | +[@roden_gedney2000] uses the complex frequency-shifted (CFS) |
| 125 | +stretching function and implements the PML via recursive |
| 126 | +convolution. The key update for the auxiliary $\Psi$ fields is: |
| 127 | +$$ |
| 128 | +\Psi_x^{n+1} = b_x \Psi_x^n + a_x \frac{\partial E_z}{\partial x}\bigg|^n, |
| 129 | +$$ {#eq-em-cpml-psi} |
| 130 | +where the coefficients $b_x$ and $a_x$ are derived from the CFS-PML |
| 131 | +parameters: |
| 132 | +$$ |
| 133 | +b_x = e^{-(\sigma_x/\kappa_x + \alpha_x)\Delta t}, \quad |
| 134 | +a_x = \frac{\sigma_x}{\kappa_x(\sigma_x + \kappa_x \alpha_x)}(b_x - 1). |
| 135 | +$$ |
| 136 | +
|
| 137 | +The CPML modifies the spatial derivatives in both the H and E updates: |
| 138 | +$$ |
| 139 | +\left.\frac{\partial E_z}{\partial x}\right|_{\text{PML}} |
| 140 | += \frac{1}{\kappa_x}\frac{\partial E_z}{\partial x} + \Psi_x. |
| 141 | +$$ |
| 142 | +
|
| 143 | +This approach is superior to the simple conductivity method because: |
| 144 | +
|
| 145 | +- It provides true impedance matching at the PML interface |
| 146 | +- It handles evanescent waves (with $\alpha > 0$) |
| 147 | +- It is stable for long-time simulations and dispersive media |
| 148 | +
|
| 149 | +Typical parameters: |
| 150 | +
|
| 151 | +- `pml_width`: 10--20 grid cells |
| 152 | +- `order`: 3--4 (polynomial order for $\sigma$ grading) |
| 153 | +- `sigma_max`: Computed from optimal formula |
| 154 | +
|
| 155 | +The optimal $\sigma_{\max}$ minimizes total reflection (numerical + PML): |
| 156 | +$$ |
| 157 | +\sigma_{\max} = \frac{(m+1)}{150\pi \Delta x} |
| 158 | +$$ {#eq-em-pml-sigma} |
| 159 | +
|
| 160 | +where $m$ is the polynomial order. |
| 161 | +
|
| 162 | +### Visualizing 2D Propagation |
| 163 | +
|
| 164 | +```python |
| 165 | +import matplotlib.pyplot as plt |
| 166 | +from src.em import solve_maxwell_2d, gaussian_source_2d |
| 167 | +import numpy as np |
| 168 | +
|
| 169 | +# Point source excitation |
| 170 | +def source(t): |
| 171 | + from src.em import ricker_wavelet |
| 172 | + return ricker_wavelet(np.array([t]), f0=1e9)[0] |
| 173 | +
|
| 174 | +result = solve_maxwell_2d( |
| 175 | + Lx=0.5, Ly=0.5, Nx=200, Ny=200, T=2e-9, CFL=0.5, |
| 176 | + source_func=source, |
| 177 | + source_position=(0.25, 0.25), |
| 178 | + pml_width=20, |
| 179 | + save_history=True, |
| 180 | + save_every=10, |
| 181 | +) |
| 182 | +
|
| 183 | +# Plot snapshots |
| 184 | +fig, axes = plt.subplots(2, 3, figsize=(12, 8)) |
| 185 | +for ax, i in zip(axes.flat, range(0, len(result.E_history), len(result.E_history)//6)): |
| 186 | + im = ax.imshow(result.E_history[i].T, origin='lower', |
| 187 | + extent=[0, 0.5, 0, 0.5], cmap='RdBu', |
| 188 | + vmin=-1, vmax=1) |
| 189 | + ax.set_title(f't = {result.t_history[i]*1e9:.2f} ns') |
| 190 | + ax.set_xlabel('x [m]') |
| 191 | + ax.set_ylabel('y [m]') |
| 192 | +plt.tight_layout() |
| 193 | +``` |
| 194 | +
|
| 195 | +The simulation shows: |
| 196 | +
|
| 197 | +1. **Circular wavefront**: Expanding from the point source |
| 198 | +2. **No visible reflections**: PML absorbs outgoing waves |
| 199 | +3. **Correct propagation speed**: Wavefront radius = $c \times t$ |
| 200 | +
|
| 201 | +### Scattering from Objects |
| 202 | +
|
| 203 | +A classic FDTD application is computing scattering from objects. We can |
| 204 | +embed a dielectric or conducting scatterer: |
| 205 | +
|
| 206 | +```python |
| 207 | +from src.em.materials import create_cylinder_model_2d, GLASS |
| 208 | +
|
| 209 | +# Create cylinder scatterer |
| 210 | +eps_r, sigma = create_cylinder_model_2d( |
| 211 | + Nx=200, Ny=200, Lx=0.5, Ly=0.5, |
| 212 | + center=(0.25, 0.25), |
| 213 | + radius=0.03, |
| 214 | + cylinder_material=GLASS, |
| 215 | +) |
| 216 | +
|
| 217 | +result = solve_maxwell_2d( |
| 218 | + Lx=0.5, Ly=0.5, Nx=200, Ny=200, T=2e-9, |
| 219 | + eps_r=eps_r, |
| 220 | + source_func=source, |
| 221 | + source_position=(0.1, 0.25), # Source to left of cylinder |
| 222 | + pml_width=20, |
| 223 | + save_history=True, |
| 224 | +) |
| 225 | +``` |
| 226 | +
|
| 227 | +The scattered field shows: |
| 228 | +
|
| 229 | +- **Reflection** from the front surface |
| 230 | +- **Transmission** through the cylinder |
| 231 | +- **Diffraction** around the edges |
| 232 | +- **Internal resonances** for certain sizes |
| 233 | +
|
| 234 | +### Connection to Wave Chapter ABCs |
| 235 | +
|
| 236 | +The PML can be viewed as a sophisticated extension of the absorbing boundary |
| 237 | +conditions discussed for the scalar wave equation (@sec-wave-abc). Compare: |
| 238 | +
|
| 239 | +| Method | Principle | Angle Dependence | Complexity | |
| 240 | +|--------|-----------|------------------|------------| |
| 241 | +| First-order ABC | One-way wave equation | Normal incidence only | Simple | |
| 242 | +| Higher-order ABC | Multiple angles | Improved | Moderate | |
| 243 | +| **PML** | Impedance matching | All angles | Higher | |
| 244 | +
|
| 245 | +The PML achieves angle-independent absorption through the impedance matching |
| 246 | +condition, which ensures zero reflection at the PML interface regardless of |
| 247 | +incidence angle. |
| 248 | +
|
| 249 | +### Practical Considerations |
| 250 | +
|
| 251 | +**Grid Resolution**: |
| 252 | +
|
| 253 | +The rule of thumb is 10-20 grid points per wavelength for acceptable accuracy. |
| 254 | +For a 1 GHz signal in vacuum: |
| 255 | +$$ |
| 256 | +\lambda = \frac{c}{f} = \frac{3 \times 10^8}{10^9} = 0.3 \text{ m} |
| 257 | +$$ |
| 258 | +
|
| 259 | +So we need $\Delta x \leq 0.03$ m (30 mm) for 10 points per wavelength. |
| 260 | +
|
| 261 | +**Computational Cost**: |
| 262 | +
|
| 263 | +2D FDTD scales as $O(N_x \times N_y \times N_t)$. For a $200 \times 200$ grid |
| 264 | +with 1000 time steps, this is 40 million field updates. Each requires only |
| 265 | +a few floating-point operations, making FDTD very efficient. |
| 266 | +
|
| 267 | +**Memory**: |
| 268 | +
|
| 269 | +We need to store 3 field arrays ($E_x$, $E_y$, $H_z$) plus material property |
| 270 | +arrays. For a $200 \times 200$ grid in double precision: |
| 271 | +$$ |
| 272 | +\text{Memory} \approx 6 \times 200 \times 200 \times 8 \text{ bytes} \approx 2 \text{ MB} |
| 273 | +$$ |
| 274 | +
|
| 275 | +This is modest by modern standards, allowing much larger simulations. |
0 commit comments