简介：面向电磁场、声学与热传导等波动问题的 Python 与 CUDA 加速 FDTD（时域有限差分）学习项目，适合具备一定 Python 基础、希望理解数值求解偏微分方程并借助 GPU 提速的读者。压缩包共 35 个文件，以 Python 脚本为主（20 个 py），另有 2 个 CUDA 内核文件、9 个 txt 物理量与说明文件，以及实现文档和开发历程记录；整体仅 1.19MB，结构紧凑，已有 183 人学习浏览。内容覆盖一维到二维的多种仿真示例，如波导、分束器、环形振荡器、PML 吸收边界、Ricker 子波源等，并包含 dipole 的 CUDA 加速实现；按一维变化、二维结构、广义坐标、正交柱坐标等主题分类，便于按需查阅。配套文本给出 h-mu、h-phi 等物理量与坐标定义，可对照公式理解 FDTD 的离散化、时间步进、边界处理，以及 GPU 线程分配、内存管理、并行效率等优化要点。下载后可快速复现仿真场景，也能在此基础上扩展自己的电磁场或声学模拟任务，适合作为学习 FDTD 与 PyCUDA 的参考素材。

1. 一个“下载.zip”背后的完整电磁仿真方案：Python + CUDA 时域有限差分法

拿到带“Python_Cuda_下载.zip”后缀的项目，懂行的人第一反应不是急着解压，而是先问：这套时域有限差分法（FDTD）代码包，给到我手里的到底是一个能直接跑的求解器，还是一堆需要自己拼装的零件？FDTD 是电磁场数值模拟里最成熟的时域算法之一，它把麦克斯韦方程组在空间和时间上交错离散，直接模拟波的传播、散射与耦合过程。纯 Python 写一个一维教学版 FDTD 不难，但网格一旦铺到百万级，解释器逐层循环会慢到让人怀疑人生，CUDA 加速正好卡在这个痛点上。这篇笔记面向准备拿这套代码做光学、天线或射频仿真的工程师和学生，从原理、最小实现到 CUDA 迁移路径，把每一步的选型理由和常见坑位一次讲透。

2. 时域有限差分法先立住原理：Yee 网格、交替更新与 Courant 条件

2.1 为什么 FDTD 能在电磁仿真里当“万金油”

FDTD 能成为工程界的常青树，靠的不是精度天下第一，而是它在时域正面求解的独特优势。频域方法通常一次求解只给一个频率点的响应，想拿到宽频结果就要扫频，计算量线性叠加；FDTD 则是对一个宽带脉冲源做时间步进，一次算完就能通过傅里叶变换得到整个频段的响应。做天线、滤波器、超材料这类强色散、强谐振问题，这个特性省下的时间不是一倍两倍。

另一个被低估的优点是对非线性介质和时变材料的天然适配。很多频域方法的核心假设是线性时不变，遇到光学非线性效应只能绕路；FDTD 的迭代结构里，每个网格点的介电常数可以随局部场强变化当场更新，非线性效应只是多写一行赋值的事。加上 FDTD 处理复杂几何只需要网格填充，不需要像有限元那样生成高质量网格，程序结构的鲁棒性也高出一截。

代价同样明显：网格分辨率直接决定了内存和算力。工程经验是每个波长至少铺 10 到 20 个网格，三维场景下网格数是波长分辨率的立方增长。一个 10 微米见方的近红外光子器件，在 1550 nm 波长下粗略也要铺 100 乘 100 乘 100 的网格，如果时间步再被 Courant 条件卡住，步进上万步是常态。这就是 CUDA 进来的理由，搞不定性能，FDTD 只能停留在教学玩具阶段。

2.2 Yee 晶格与 E/H 场的交错离散：代码里那两行更新公式从哪来

FDTD 最核心的设计出自 Kane Yee 在 1966 年的论文：把电场和磁场在空间上错开半个网格，在时间上错开半个步长。以三维为例，每个电场分量周围环绕四个磁场分量，每个磁场分量周围同样环绕四个电场分量。这个交错布局让麦克斯韦旋度方程可以被二阶中心差分直接离散，不用解任何线性方程组，当前时刻的场值只由上一时刻的相邻场值决定。

一维情况下，假设波沿 x 方向传播，只保留 Ez 和 Hy 两个分量，离散位置写成：

• Ez 定义在整数网格点 i 上
• Hy 定义在半整数网格点 i + 0.5 上

从法拉第定律和安培定律出发，展开旋度方程后得到：

# 安培定律：E 的更新用 H 在相邻两个半网格点的差值 ez[i] = ez[i] + courant * (hy[i] - hy[i-1]) # 法拉第定律：H 的更新用 E 在相邻两个整数网格点的差值 hy[i] = hy[i] + courant * (ez[i+1] - ez[i])

courant 在这里是库朗数，真空归一化单位下等于 c 乘以 dt 除以 dx。这套交错更新的意义在于，E 场的旋度只依赖 H 场的空间差分，H 场同理，每个网格点上的值独立更新，天然适合并行。三维场景只是重复六次相同模式的更新，分别对应 Ex、Ey、Ez 和 Hx、Hy、Hz，索引错位的套路完全一致。

真正容易写错的是索引边界。Ez 的更新用到了 hy[i] 和 hy[i-1]，那么 Ez 最左边一个格子依赖不存在的 hy[-1]；Hy 的更新用到了 ez[i+1]，最右边一个格子同样越界。代码里必须把两端单独处理，或者把更新的索引范围收窄，不然数组越界在 Python 里可能不报错，直接用垃圾值污染所有结果，这也是新手最容易翻车的地方。

2.3 Courant 稳定性条件与网格分辨率：参数设错的典型翻车现场

FDTD 是显式时间步进算法，时间步长不能随便取。稳定性条件由库朗数约束：三维情况下 courant 必须小于等于 1 除以根号 3，二维是 1 除以根号 2，一维就是小于 1。实际使用时通常取 0.9 倍的安全系数，三维用 0.5 到 0.7 都是常见做法。这个限制来源于差分格式的特征值谱，超了之后误差会随时间步指数放大，表现出来就是波形跑到几百步突然变成 NaN。

比发散更隐蔽的是数值色散。即使满足稳定性条件，时间和空间差分近似本身会引入与频率相关的相速度误差。网格分辨率越低，误差越大。对宽频脉冲而言，不同频率成分的波会以略有差异的速度传播，导致波形在传播几个波长后展宽变形。这种数值色散是个“玄学”问题，不报错但污染结果，唯一可靠的对策就是把网格取到每波长 20 个网格以上，然后在关键算例里做分辨率收敛性验证。

网格分辨率还直接和内存挂钩，这里给一个预估值。一个 256 立方网格的单精度场数组，6 个分量（Ex、Ey、Ez、Hx、Hy、Hz）大约占用 400 MB，双精度翻倍。如果网格翻到 512 立方，单精度也要 3.2 GB。加上 PML 吸收边界和源区网格，趁早用单精度能省一半显存，这在高分辨率三维模拟里往往是成败手。

3. 用 Python 跑通一维 FDTD：最小代码、参数含义与结果自检

3.1 代码骨架：源项注入、场更新、边界处理三步走

一套可复现的 FDTD 主循环，不管几维，都绕不开四个环节：更新 E 场、注入源项、更新 H 场、吸收边界。顺序不能乱，因为 E 和 H 在时间上差半步，先算谁后算谁决定了时间对齐关系。常规做法是先更新 E（用上一时刻的 H），然后注入源，再更新 H（用刚算好的 E），这样源项对 E 的影响会在下一步传递到 H，物理上是对的。

源项的类型决定了模拟的用途。高斯脉冲是默认选择，因为它的频谱是平滑的高斯形，覆盖从低频到高频的连续范围，适合宽频响应分析；正弦连续波适合单频稳态问题，比如谐振腔的 Q 值计算；调制高斯脉冲则兼顾两者，适合需要频带限制的场景。做散射问题还可以用平面波源加总场散射场分区，不过那属于进阶玩法，先把点源跑通再升级。

边界处理是最容易被忽略的一部分，初学者经常用零边界或者硬边界，结果波一到边界就弹回来，在域内来回反射，模拟结果被污染得面目全非。我一般建议从一阶 Mur 边界起步，它只需要一个赋值语句把边界格子的场值设为相邻格子场值，本质上是行波假设下的简易无反射截断，对教学和小规模验证足够。追求低反射率就用 PML，这留到最后一章展开。

3.2 一维自由空间传播的完整脚本与逐段说明

import numpy as np # 网格与物理参数 nx = 800 # 空间网格数 dx = 1e-3 # 空间步长，单位米 courant = 0.8 # 库朗数，必须小于 1 dt = courant * dx / 3e8 # 时间步长，单位秒 nt = 600 # 总步进数 # 场初始化 ez = np.zeros(nx) hy = np.zeros(nx) # 高斯脉冲源参数 src_pos = nx // 2 # 源放在域中间 t0 = 80 # 脉冲峰值时刻 spread = 20 # 脉冲宽度控制 # 主循环 for n in range(nt): # 第一步：安培定律更新 Ez for i in range(1, nx): ez[i] = ez[i] + courant * (hy[i] - hy[i-1]) # 第二步：注入高斯脉冲源 ez[src_pos] += np.exp(-((n - t0) / spread) ** 2) # 第三步：法拉第定律更新 Hy for i in range(nx - 1): hy[i] = hy[i] + courant * (ez[i+1] - ez[i]) # 第四步：一阶 Mur 吸收边界 ez[0] = ez[1] ez[-1] = ez[-2] # 每 100 步打印一次峰值，确认场没有发散 if n % 100 == 0: print(f"step={n}, max|Ez|={np.max(np.abs(ez)):.4f}")

这段代码的核心逻辑就是上一章的两条更新公式套进循环。courant 取 0.8 是因为一维稳定性上限是 1.0，留 20% 余量既保证数值稳定，又不会因为步长太小拖慢完成时间。dt 由 courant 除以光速换算而来，这一步不能省，后续做有耗介质或色散介质时会遇到物理时间相关信息，提前把这个换算关系写清楚省得到时候对着步数换算物理时间。

源项注入放在 E 更新之后、H 更新之前，时间对齐顺序前面已经说明。高斯源的表达式里 t0 控制峰值出现的步数，spread 控制脉冲时域宽度，他们共同决定了激励信号的频谱覆盖范围。想激励高频分量就缩小 spread，但网格分辨率跟不上时高频分量会被数值色散歪曲，这是一个需要协作调节的参数对。

Mur 边界的实现就是把边界格子的场值复制为内侧相邻格子的值。这个操作相当于假装波从边界“走出去”，对垂直入射的波有效，斜入射时反射系数会明显上升。这也是为什么后面要做 PML 的根本原因，一阶 Mur 应付教学场景可以，做工程算例必须换方案。

3.3 怎么确认代码没写错：能量守恒与波形形态两个判据

代码能跑通是一回事，跑出的结果是对的是另一回事。FDTD 代码写好后，第一步不是上大网格，而是先用这个一维脚本做两件事：观察波形是否按预期向两侧传播，以及检查能量是否守恒。两个判据都能在代码里以很低的成本实现。

能量守恒判据的做法是每步计算电磁场总能量，公式是电场能量密度加磁场能量密度的空间积分。在一维自由空间里，没有源和损耗时总能量应该基本不变，只会在源注入期间上升，随后保持平稳。如果能量一直在涨，说明 Courant 条件被破坏；如果能量突然掉头，检查边界和网格是否异常。

# 在每个时间步末尾累加，可以发现早期发散 total_energy = 0.5 * np.sum(ez**2) + 0.5 * np.sum(hy**2)

波形形态的判据更好操作：运行模拟后取某一时刻的 ez 快照，应该看到源点两侧各有一个对称的脉冲向外传播。脉冲峰值位置大致等于源位置加减光速乘以物理时间，可以用这个位置偏差来验证 dt 换算对不对。我见过不少 dt 错算成 1 导致波形直接飞出域外的案例，这种一眼就能看穿的低级错误，跑一次波形快照比查半天代码快得多。

如果这两项都过关，这个一维骨架就可以放心作为后续功能开发的基座。加介质层、加 Mur 边界升级、加金属板反射，都是在这个循环里插代码的事。等三维版本跑不动了，再进入 CUDA 迁移阶段。

4. 从 Python 到 CUDA：三种加速路径与一条实践路线

4.1 纯 Python 循环慢在哪：三维网格就是显存和算力的黑洞

一维 FDTD 用 Python 双循环完全没问题，但三维场景下问题就来了。一个常见测试是 256 立方网格跑 1000 步，光场数组就有 6 个分量约 400 MB，纯 Python 循环每步要更新约 1600 万网格点乘 6 个分量，接近一亿次赋值操作，Python 解释器逐行执行这个量级的循环，一个算例跑到天荒地老。

用 NumPy 切片能救一部分，把内层循环换成数组运算后，单步耗时能降到毫秒级。但三维数组切片加广播引入的临时数组和内存拷贝，在显存有限的场景里依然是瓶颈。更麻烦的是 FDTD 的时间步进是串行的，每步依赖上一步完成，没法把时间循环也并行掉。所以真正的出路是把空间网格的更新并行化，而这正是 CUDA 最擅长的形态：一大块同构数据、大量线程、少量逻辑判断。

加速还得先搞清楚天花板在哪。FDTD 每个网格点的运算量很小，主要耗时来自从显存读入两个相邻场值再写回一次，属于典型的带宽密集型任务。实测经验是，一个现代 GPU 的显存带宽能达到 1 TB/s 级别，按每个网格点更新需要读 16 字节写 4 字节来算，理论吞吐量也能到百亿级网格更新每秒。相比 CPU 单核的个位数亿次，数量级差在这里。

4.2 路径一：Numba 的 cuda 目标，改动最小、上手最快

Numba 是 Python 加速里门槛最低的选项，不需要离开 Python 环境，也不需要手写 C/C++。把之前的一维循环搬进 GPU，只需要把更新函数换成@cuda.jit装饰器，再显式配置线程块和网格。代价是你必须把数组显式传进 GPU 显存，并且接受 Numba 的 CUDA 支持版本绑定的约束。

from numba import cuda import numpy as np nx = 800 nt = 600 courant = 0.8 # 在 GPU 显存上分配场数组 ez = cuda.to_device(np.zeros(nx, dtype=np.float32)) hy = cuda.to_device(np.zeros(nx, dtype=np.float32)) @cuda.jit def update_e(ez, hy, courant): i = cuda.grid(1) if 0 < i < ez.shape[0] - 1: ez[i] = ez[i] + courant * (hy[i] - hy[i-1]) @cuda.jit def update_h(ez, hy, courant): i = cuda.grid(1) if 0 <= i < hy.shape[0] - 1: hy[i] = hy[i] + courant * (ez[i+1] - ez[i]) threads_per_block = 128 blocks_per_grid = (nx + threads_per_block - 1) // threads_per_block src_pos = nx // 2 t0 = 80 spread = 20 for n in range(nt): update_e[blocks_per_grid, threads_per_block](ez, hy, courant) # 源项注入只能从 host 端访问到设备数组，先用 copy 回再写 if n % 50 == 0: # 示例：低频次注入，模拟源 host_ez = ez.copy_to_host() host_ez[src_pos] += np.exp(-((n - t0) / spread) ** 2) ez = cuda.to_device(host_ez) update_h[blocks_per_grid, threads_per_block](ez, hy, courant)

Numba 分两个 kernel 更新 E 和 H，而不是在同一个 kernel 内做完，是为了避开线程间同步问题。FDTD 的 E 更新依赖上一时刻 H，H 更新依赖新 E，如果放进同一个 kernel，必须在中间插入cuda.syncthreads()，而 syncthreads 要求所有线程都到达同步点，线程块内的逻辑分支不一致时会卡死。拆成两个 kernel 调用后，kernel 之间天然同步，安全省心。

一个明显的坑是源项注入。Numba CUDA kernel 内部做ez[src_pos] +=这种单点操作效率不高，而且需要把数组拷回 host 再写再传回，虽然一维场景无所谓，但三维里频繁搬运会成为性能杀手。常见做法是把源项改成点源数组，在 kernel 内部按时间步判断全局线程 ID 是否命中源位置，把注入逻辑并进 update_e 一起算。

4.3 路径二：PyCUDA 手写 kernel，控制最细、调试最痛苦

PyCUDA 把 CUDA 的 C 风格 kernel 直接塞进 Python 源码里，用运行时编译生成 GPU 机器码。它能精确控制每个线程做什么、用多少寄存器、怎么对齐显存访问，是三维大规模算例里性能最稳定的路径。代价是你要同时维护 Python 侧的数据流和 C 侧的 kernel 代码，编译报错时还得看 C 的语法。

import pycuda.autoinit import pycuda.driver as drv from pycuda.compiler import SourceModule import pycuda.gpuarray as ga import numpy as np mod = SourceModule(""" __global__ void update_e(float *ez, const float *hy, float courant, int nx) { int i = blockIdx.x * blockDim.x + threadIdx.x; if (i > 0 && i < nx - 1) { ez[i] = ez[i] + courant * (hy[i] - hy[i-1]); } } __global__ void update_h(float *ez, float *hy, float courant, int nx) { int i = blockIdx.x * blockDim.x + threadIdx.x; if (i >= 0 && i < nx - 1) { hy[i] = hy[i] + courant * (ez[i+1] - ez[i]); } } """) update_e = mod.get_function("update_e") update_h = mod.get_function("update_h") nx = 800 ez = ga.zeros(nx, dtype=np.float32) hy = ga.zeros(nx, dtype=np.float32) block = (128, 1, 1) grid = ((nx + 127) // 128, 1) for n in range(nt): update_e(ez, hy, np.float32(courant), np.int32(nx), block=block, grid=grid) # 源注入直接操作 GPU 数组 ez[src_pos] = ez[src_pos] + np.exp(-((n - t0) / spread) ** 2) update_h(ez, hy, np.float32(courant), np.int32(nx), block=block, grid=grid)

PyCUDA 的get_function返回的 kernel 对象支持block和grid配置，参数在调用时通过 Python 对象直接传入，kernel 内部按 C 的指针语义访问显存。这里比 Numba 多了一个能力：源项注入可以直接在 host 端操作 GPU 数组的一个元素，PyCUDA 会自动处理单元素读写，省去了 Numba 里的来回拷贝。这是 PyCUDA 在工程效率上的实在优势。

索引和步长的匹配是 PyCUDA 新手最常犯的错。C kernel 里的nx参数和 Python 侧的数组长度必须一致，block 大小要和网格划分配合。例子里的grid = ((nx + 127) // 128, 1)是把总网格数按每块 128 线程切分，如果 nx 不是 128 的整数倍，最后一个块会有部分线程不满足边界条件退出，这是正确的写法。很多人图省事直接grid=(10,1)导致线程数超过数组长度，越界访问是硬隔离的，直接崩掉。

4.4 路径三：C++/CUDA 编译扩展，留给性能强迫症

Numba 和 PyCUDA 的运行时编译都有小幅启动开销，kernel 本身也受限于各自代码生成器的优化能力。当网格规模上到千万级、单次仿真要跑几万步时，用 C++/CUDA 写 kernel 编译成 Python 扩展模块是性能上限最高的路线。BDF 的维护成本高，但数据结构完全自己掌控，还能用 shared memory 做相邻场值的缓存优化。

// fdtd_kernel.cu 核心更新（一维示例，三维同理扩展到三个方向） __global__ void update_e(float *ez, const float *hy, float courant, int nx) { int i = blockIdx.x * blockDim.x + threadIdx.x; if (i > 0 && i < nx - 1) { ez[i] = ez[i] + courant * (hy[i] - hy[i-1]); } }

编译命令大致长这样，注意-arch参数要匹配目标 GPU 的计算能力。以 RTX 4060 Ti 为例，它的计算能力是 8.9，对应sm_89，编译时指定这个参数让编译器生成适配指令集的机器码。用旧的默认架构编译也能跑，但性能可能差 20% 到 30%，属于白扔的算力。

nvcc -O3 -arch=sm_89 -Xcompiler -fPIC -shared fdtd_kernel.cu -o fdtd_kernel.so

拿到.so动态库后，在 Python 里用ctypes或cffi加载调用。这要求你对 C ABI 和内存生命周期有足够把握，因为数组指针要显式从 NumPy 转换并保持引用，否则 Python 侧可能把底层缓冲释放，kernel 访问到悬空指针，表现为偶发的 Segmentation Fault，这类问题极难排查。所以我的建议是：没有明确的性能瓶颈，不要上编译扩展；Numba 和 PyCUDA 对绝大多数情况已经够用。

4.5 加速比不是玄学：显存带宽、网格规模与 kernel 合并访问

很多人拿一维 FDTD 跑完 GPU 版本后抱怨加速比只有几倍，这不是 CUDA 不行，而是问题形态决定了它跑不满 GPU。一维网格只有 800 个点，一个 kernel 只启动几个线程块，GPU 的数千个计算核心大部分在空转。CUDA 的优势在大规模并行，一维场景连显存带宽都吃不满，自然拿不到好看的数字。

真正让 GPU 发挥实力的是二维和三维。二维 512 乘 512 网格单次更新就是 26 万线程，三维一个面就有百万线程，这时候并行度才够。实测数据可以参考：同样 256 立方网格跑 1000 步，纯 Python 需要小时级，NumPy 向量化大约几十分钟，PyCUDA 能压到一两分钟，数量级差距通常在 20 到 50 倍，具体数字取决于显卡型号和参数设置。

如果加速比达不到预期，优先检查内存访问合并。FDTD 更新时每个线程读相邻两个网格的值，如果线程 ID 与网格索引一一对应，相邻线程读的是相邻内存地址，这条访问路径是合并的，带宽利用率接近理论峰值。如果索引映射走了形（比如线程 ID 乘了 2），同一个 warp 的线程读的内存跨度翻倍，带宽利用率直接腰斩。另一个常见问题是数据类型，double 的访问带宽只有 float 的一半，能上 float 就别用 double，除非双精度需求硬性存在。

# 检查当前显卡计算能力，确认 kernel 能覆盖目标架构 import pycuda.driver as drv dev = drv.Device(0) print(f"compute capability: {dev.compute_capability()}")

5. 实战避坑：CUDA 环境、显存边界与数值发散的五个典型现场

5.1 现象一：CUDA out of memory，三维网格把显存吃光了

现象：代码跑得好好的，网格规模一放大，启动 kernel 时直接报CUDA out of memory，或者程序在数组分配阶段就挂掉。常见于从二维测试切到三维正式算例的时候。

原因：三维场数组的内存占用按立方增长。一个 N 立方网格，6 个场分量，float32 单精度的占用是 N 的立方乘以 24 字节；N=256 时约 400 MB，N=512 时约 3.2 GB，加上 PML 区域和源数组，8 GB 显存的卡很快就见底。这是硬件边界的数学现实，不是代码 bug。

解决：先算后跑，别开盲盒。把目标网格规模代入/ 1e9算 GB 数，低于显存一半才开工。然后优先把 double 降到 float32，很多模拟场景对精度的需求并没有高到非双精度不可；再用网格渐变或者子区域拆分把大问题化小。还有一个实用招数：只用一个场数组做滚动更新，不要为了代码美观为每个时间步保留历史快照，历史数据用磁盘流式存储。

5.2 现象二：驱动版本与 CUDA Toolkit 对不上，设备初始化直接失败

现象：pycuda.driver.Error: CUDA driver version is insufficient for CUDA runtime version，或者 Numba 报cuda.cuda_error.CudaSupportError。看起来像代码问题，实际是环境不匹配。

原因：NVIDIA 驱动向后兼容，新驱动能跑旧版 CUDA runtime，但旧驱动跑不了新版 runtime。比如驱动 550.144.03 支持 CUDA 12.x，如果你装了 CUDA 12.6，没问题；但环境里残存着 CUDA 13 的 runtime 库，就必然报 insufficient。这是加速器生态环境里最常见的“装完就跑不通”场景。

解决：先看硬件支持的 CUDA 上限，再看代码依赖的上限，最后取交集。命令行直接查：

nvidia-smi # 看右上角 CUDA Version，这是驱动支持的最大 runtime 版本 nvcc --version # 看当前激活的 CUDA Toolkit 版本

Toolkit 版本必须小于等于驱动支持版本。如果不得已要多版本共存，用conda install cudatoolkit装指定版本的 CUDA 运行库，和系统的nvcc解耦，能解决九成环境冲突。Numba 用户还要额外检查numba.__version__对应的 CUDA 绑定，Numba 0.58 绑 CUDA 11.8，Numba 0.59 绑 CUDA 12.x，版本错位会出现找不到运行时函数的诡异报错。

5.3 现象三：WSL2 里找不到设备，CUDA 在 Linux 子系统里的特殊要求

现象：在 WSL2 的 Ubuntu 里 import pycuda 或 numba.cuda，报错说找不到 CUDA 设备；nvidia-smi命令也不存在，但 Windows 侧驱动工作正常。

原因：WSL2 的 CUDA 架构和原生 Linux 不一样。它复用 Windows 侧安装的 NVIDIA 驱动，WSL 内部不需要也不能直接安装单独的驱动内核模块。很多人按 Linux 教程在 WSL 里装了一堆驱动相关包，装完反而把驱动版本搞乱，设备还是起不来。

解决：按 WSL 专用路径走。第一步在 Windows 侧安装 NVIDIA WSL 驱动（不是标准桌面版驱动，是带 WSL 支持的版本），第二步在 WSL 内部只装 CUDA Toolkit 用户态部分，不要碰内核驱动。验证方式是：

nvidia-smi # 能看到 GPU 型号和驱动版本，但 CUDA 版本可能显示 12.x python -c "import pycuda.autoinit; print(pycuda.driver.Device(0).name())"

WSL2 下跑 CUDA 的性能接近原生 Linux，但 PCIe 直通和显存访问的延迟会比原生高一些，对 FDTD 这种带宽敏感型任务影响有限。如果跑出来性能明显不对劲，优先检查是不是误用了 Windows 侧的路径和库。

5.4 现象四：波形发散成 NaN，Courant 条件被破坏

现象：模拟跑到几百步后，场值逐渐变大然后变成 NaN，或者直接一步跳变。打印每个时间步的能量曲线，能看到指数增长。

原因：最常见的是 courant 大于 1。三维情况下稳定性上限是 1 除以根号 3，约 0.577，很多人直接沿用一维 0.8 的经验值，直接爆了。另一个隐蔽原因是网格里包含不同介质区域，介质的折射率不同导致局部光速不同，如果 courant 是在真空参数下算的，在折射率大于 1 的介质里 dt 不变但局部 c 变小，等效库朗数下降，这不会导致不稳定但会让误差变大。反过来在折射率小于 1 的色散介质里就可能突破稳定上限。

解决：全局按最坏情况取 courant。三维固定取 0.5，二维取 0.7，一维取 0.8，这是个不容易出错的规则。然后每 10 步打印一次np.max(np.abs(ez))和总能量，看到数值开始单调上涨就立刻停，用二分法找临界步长，比一次跑完看结果再猜高效得多。

还有个被人忽视的细节：边界条件写成了硬边界而不是吸收边界时，波会在域内来回反射并叠加，虽然数值本身不发散，但峰值会明显虚高，波形永远平稳不下来。这种看起来不像报错的“隐性翻车”，只能靠能量守恒判据来抓。

5.5 现象五：边界反射把结果污染了，缺了吸收边界

现象：波形在传播到模拟区域边界后反射回来，在域内形成驻波状干涉条纹，反射波的峰值甚至能和原始波形叠加出错误的共振峰。

原因：一阶 Mur 边界只对垂直入射波有效，斜入射反射系数随角度增大，三维模拟中波从任意角度打到边界的情况太普遍，Mur 边界撑不住。更关键的是如果干脆把边界设成零值硬边界（PEC 默认行为），波会在边界全反射，数值上有解但物理上完全不是你想模拟的场景。

解决：上 PML。PML 的思路是在边界区域植入一层有损耗介质，让进入的波指数衰减，到最外层再反射回来时已经衰减到可以忽略的程度。一个设计良好的 PML，理论反射系数能到 1e-6 以下，远好于 Mur 边界的百分之一量级。代价是边界区域要多算十几个网格层的更新，并引入额外的辅助变量，内存和算力多消耗约 20% 到 30%。

我测试过最直接的反差：同样的 100 乘 100 二维 FDTD，用 Mur 边界跑出来的透射率曲线在高频段有明显波纹，换成 10 层 PML 后波纹消失，曲线干净得像教科书插图。这个改善不值得犹豫，从工程第一天就应该用 PML。

6. 把这个代码包用出价值：PML 吸收边界与近远场外推的落地技巧

PML 是这个代码包里最值得先吃透的模块，因为几乎所有真实算例都依赖它。它的数学形式看起来复杂，实现上却有一层窗户纸。以 x 方向 PML 为例，要点是在 PBL 区域内把每个电场分量的更新公式改成带损耗系数的形式，损耗系数 sigma 从边界内沿向最外层按多项式梯度递增。常见做法是 10 层厚度、4 次方坡度，理论反射率在 1e-5 以下。

# 一维 PML 区域中 Ez 的更新，sig 按层数梯度分布 # ca 和 cb 是最重要的参数，拿回来和总是 1 的自由空间公式对照 sig = sigma_max * (x_in_pml / pml_thickness) ** 4 ca = (1 - dt * sig / 2) / (1 + dt * sig / 2) cb = courant / (1 + dt * sig / 2) ez[i] = ca[i] * ez[i] + cb[i] * (hy[i] - hy[i-1])

ca 项逐渐把旧场值压下去，cb 项仍然传递相邻网格的耦合，两者配合让波在层层推进中不断衰减。调参的核心是 sigma_max：取太小衰减不够，取太大数值反射反而上升，工程上先在自由空间场景里测一遍反射率，用峰值反射系数作为调参目标，再投入正式算例。

近远场外推是让这个代码包从“能算场分布”升级到“能算远场方向图”的关键一步。FDTD 天然算的是近场区域的场值，而天线远场方向图、雷达散射截面这些工程指标是在远场定义。常规做法是：在包围目标的闭合曲面上记录等效电流和等效磁流，用等效原理在频域做积分变换。这个过程不复杂，但要记得在时间步进结束后把所有记录场值做一次傅里叶变换，再到关心的频率点上执行外推计算。

拿这套代码做工程项目，我最想提醒的是控制预期。别一上来就铺三维大网格加 PML 加近远场外推，这套链路里每一环都有一到两个隐藏参数，任何一个设错都让结果失真。最好的习惯是先跑通一维透射验证，再上二维环形谐振腔，最后才碰三维天线模型。给自己写一个回归脚本：模拟自由空间高斯源传播，断言波形峰值位置与理论位置偏差不超过一个网格，改任何代码后先跑这个脚本，全绿之前不碰新功能。这个习惯帮我躲过了无数个从下午修到凌晨的 debug 局。希望帮到你。

本文还有配套的精品资源，点击获取