简介：一套基于C语言的有限差分时域法（FDTD）并行计算实现，面向计算电磁学方向的开发者与研究者，可用于模拟电磁波传播、天线辐射等场景。压缩包内共26个文件，以.h头文件和.cpp源文件为主，另有txt配置、docx文档、bat批处理等，包体约246KB，结构紧凑，便于直接阅读和编译调试。已有439人浏览学习。项目实现了点源设置、电场磁场迭代更新、边界条件处理，并通过OpenMP并行化加速大型网格计算，同时附有可执行程序与构建脚本，支持在配置文件中调整点源参数，适合入门并行FDTD算法或作为二次开发的代码基底。
1. 有限差分时域法并行C代码：先解决算得动，再解决算得准
仿真一个微带天线，网格规模一亿，单核C程序可能要跑三周。有限差分时域法（FDTD）把麦克斯韦旋度方程在时间和空间上差分离散，每个时间步只更新邻域格点，天然适合做并行；也正因为它只更新邻域，单核算力就成了瓶颈。并行化在计算电磁学里不是锦上添花，而是把仿真周期从周变成小时的常规手段。C语言在这个领域依然是主力，原因很直接：指针和数组布局能精确控制内存，MPI/OpenMP的接口就是C绑定，编译器对连续循环的自动向量化效果也最好。下文从Yee网格的核心更新方程说起，给出一套C语言的MPI并行骨架，并把CFL条件、守护单元、非阻塞通信、负载均衡这些决定成败的参数和坑位说清楚。这套路径同样适合用现代C写嵌入式场仿真或异构移植的工程师，因为并行的正确性验证方法和参数边界是一致的。
2. 有限差分时域法的Yee网格与C语言更新方程
2.1 电场磁场错开半个网格，C数组为什么按分量存放
Yee网格把电场和磁场分量错开半个网格步长，时间上也错开半个步长。这样做的结果是，电场和磁场各自用相邻点的空间差分去推进另一场，空间差分在交错点上是中心对称的，因而达到二阶精度；时间上采用蛙跳格式，也是二阶精度。如果没有这个错位，用常规网格的中央差分容易引入晶格色散，精度也会掉到一阶。C语言实现Yee网格的常见方式是按分量分数组，而不是把六个分量放进一个结构体数组。原因在缓存：FDTD一次更新只读写一两个分量，结构体数组会让一次缓存行携带六个分量，只有两个被使用，带宽浪费明显；按数组结构则每个缓存行全是同一分量，循环向量化和MPI非阻塞通信也都好处理。
2.2 磁场与电场更新循环的C代码骨架
二维TM模式虽然简单，但三维代码的循环结构完全一致，只是额外增加维度。下面这段代码展示了磁场Hx和电场Ez的核心更新：
#define IDX(i, j, ny) (((i) * (ny)) + (j)) // 磁场 Hx 更新：Hx -= (dt/(mu*dy)) * (Ez[i][j+1] - Ez[i][j]) void update_hx(double *restrict hx, const double *restrict ez, int nx, int ny, double coef_hx) { for (int i = 0; i < nx; i++) { const double *ez_row = ez + i * ny; double *hx_row = hx + i * ny; for (int j = 0; j < ny - 1; j++) { hx_row[j] -= coef_hx * (ez_row[j + 1] - ez_row[j]); } } } // 电场 Ez 更新：Ez += (dt/(eps*dx)) * (Hy[i+1][j] - Hy[i][j]) void update_ez(double *restrict ez, const double *restrict hy, int nx, int ny, double coef_ez) { for (int i = 0; i < nx - 1; i++) { const double *hy_row = hy + i * ny; const double *hy_next_row = hy + (i + 1) * ny; double *ez_row = ez + i * ny; for (int j = 0; j < ny; j++) { ez_row[j] += coef_ez * (hy_next_row[j] - hy_row[j]); } } }这里coef_hx和coef_ez是提前算好的系数，分别等于dt/(mu0*dy)和dt/(eps0*dx)，把除法移出最内层循环。IDX宏只是示意，性能敏感版本里应像上面这样用行指针递增，避免每次循环都做乘法和加法。restrict告诉编译器hx和ez没有重叠，编译器才能放心向量化。二维TM版还需要update_hy，把Hy用Ez的空间差分更新，再按同样模式扩展成三维的六个分量。三维时最内层循环仍在j方向连续，i和k放到外层，这样才能让缓存命中率最高。
2.3 吸收边界PML的参数与C数据布局
吸收边界是另一个绕不开的话题。PML（完美匹配层）通常放在计算域最外面8到10层，电导率从内到外按多项式渐变：sigma_x(x)=sigma_max*(x/d)^n，n取3或4，d是PML厚度。工程上常取理论反射率R0=1e-6，sigma_max=-(n+1)*ln(R0)/(2*eta*d)，其中eta是介质波阻抗。C语言实现时我习惯在数组四周统一加pml层，让内部网格索引都加上一个偏移量，这样更新循环内不需要判断当前是否在PML区域，只需要在时间步里把PML系数数组也一起更新。注意每个格点对应的eps和sigma都需要独立数组，电导率在PML区域内逐点变化。这些数组用一次malloc分配、按一维连续方式索引，这样后续做MPI域分解时，可以用MPI_Type_vector准确描述子域边界的非连续内存切片，通信代码才写得干净。
2.4 CFL稳定条件与并行无关，但决定总步数
网格步长和时间步长不是独立的。对于均匀介质，CFL条件要求c*dt <= 1/sqrt(1/dx^2 + 1/dy^2 + 1/dz^2)。三维均匀网格下就是dt <= h/(c*sqrt(3))。实际程序里取CFLN=0.9，即时间步长为上限的90%。如果目标频率f对应波长lambda_min，先按每波长10到20个网格选dx，再按CFL选dt。并行化不改这个条件，总步数和计算量固定，所以并行的意义是把单步的墙钟时间压下来。常用参数组合如下表。
| 参数 | 推荐区间 | 说明 |
|---|---|---|
| 网格步长dx | lambda_min/10 ~ lambda_min/20 | 决定分辨率和PML厚度 |
| CFL安全系数 | 0.8 ~ 0.9 | 过高不稳，过低增加步数 |
| PML厚度 | 8 ~ 12格 | 太少反射明显，太多浪费内存 |
| PML渐变阶数n | 3 ~ 4 | n越大，高频吸收越好 |
这个表格在并行FDTD里还有一个作用：PML区域位于全局计算域最外侧，只有部分进程需要实际分配PML电导率数组，否则每个子域都按完整PML厚度分配，内存和通信都会翻倍。
3. MPI并行化：计算域划分、守护单元与信息交换
3.1 为什么选空间域分解
FDTD的更新格式只读取相邻格点的值，计算域内任意两点相距超过一个网格时没有直接依赖，这给空间域分解提供了理论依据。将计算域沿一个或多个方向切成子域，每个MPI进程负责一块；每步迭代前，子域边界上的一层格点需要从相邻进程拿，这就是所谓的守护单元格（ghost cells）。时间域并行也有研究，但双曲型方程在时间方向的长程相关性让校正步复杂化，工程上几乎没有项目敢直接上，所以这里只讲空间分解。
3.2 用MPI_Dims_create和MPI_Cart_create构造二维虚拟拓扑
并行FDTD最常用的拓扑是笛卡尔网格。用MPI自带函数创建虚拟拓扑，比自己算邻居rank再小心处理周期条件可靠得多。
#include <mpi.h> void build_cart_comm(MPI_Comm *cart_comm, int *nbr_left, int *nbr_right, int *nbr_down, int *nbr_up) { int size; MPI_Comm_size(MPI_COMM_WORLD, &size); int dims[3] = {0, 0, 0}; int ndims = 2; // 2D解析；3D改为3 MPI_Dims_create(size, ndims, dims); int periods[3] = {0, 0, 0}; // 非周期边界 MPI_Cart_create(MPI_COMM_WORLD, ndims, dims, periods, 1, cart_comm); MPI_Cart_shift(*cart_comm, 0, 1, nbr_down, nbr_up); MPI_Cart_shift(*cart_comm, 1, 1, nbr_left, nbr_right); }MPI_Dims_create传入0会让库自动选数值，尽量让维度接近，使得子域表面最小，通信量最小。MPI_Cart_create最后一个参数是reorder，设为1允许MPI重新排列进程rank以匹配底层节点拓扑，对共享内存节点和跨节点混跑都有帮助。MPI_Cart_shift返回轴上的前后邻居；对于非周期网格，边界方向上的邻居是MPI_PROC_NULL，发送接收时要专门处理，否则MPI调用仍然成功，逻辑上却容易出错。
3.3 非阻塞通信与守护单元格交换
有了拓扑，接下来是每个时间步在子域边界交换场分量。以二维为例，每条边需要交换的数据条数等于边界上格点数。使用MPI_Isend/MPI_Irecv配合MPI_Waitall。基本通信结构如下：
// 假设 ez 需要与左右邻居交换，edge 是一条边的数据字节数 int edge = ny_local * sizeof(double); MPI_Request reqs[4]; // 先注册接收 MPI_Irecv(ez_ghost_xmin, edge, MPI_BYTE, nbr_left, tag, cart_comm, &reqs[0]); MPI_Irecv(ez_ghost_xmax, edge, MPI_BYTE, nbr_right, tag, cart_comm, &reqs[1]); // 填充发送缓冲区并发送 memcpy(send_xmin, &ez[IDX(0, 0, ny_local)], edge); memcpy(send_xmax, &ez[IDX(nx_local - 1, 0, ny_local)], edge); MPI_Isend(send_xmin, edge, MPI_BYTE, nbr_left, tag, cart_comm, &reqs[2]); MPI_Isend(send_xmax, edge, MPI_BYTE, nbr_right, tag, cart_comm, &reqs[3]); MPI_Waitall(4, reqs, MPI_STATUSES_IGNORE);这里先Irecv再Isend可以避免通信双方都阻塞在Send上的经典死锁，同时让发送缓冲区在整个Waitall之前不被复用。memcpy源地址要按子域实际布局写，ny_local是子域列数，nx_local是行数，示例中IDX(0,0,ny_local)对应子域最左边界的第一行；三维时还要把ny_local换成ny_local*nz_local，并且为每个场分量做一遍同样过程。别把ghost区域本身作为发送源，否则会把上一轮的旧值发出去。
3.4 通信量估算与检查点里的C文件操作
同样全局网格，一维分解的通信面是二维切块的若干倍。进程数较多时，二维或三维切块几乎是必选。
| 分解方式 | 二维全局规模 | 单步通信量（double） |
|---|---|---|
| 一维切条 | Nx×Ny, P行 | 2×Ny |
| 二维切块 | Nx×Ny, Px×Py | 2×(Ny/Py+Nx/Px) |
| 三维切块 | Nx×Ny×Nz | 3个面的面积之和 |
紧凑分解不仅减少通信总量，也减少每个子域的守护单元格数量。负载均衡上，每个进程的子域体积应当相同；MPI_Dims_create已经保证维度数接近，但网格维度不能被进程数整除时，各子域会差一行。常见做法是进程数拆成可整除的因子，比如64进程选2×32还是4×16，要量一下边界通信与MPI_Type_vector的代价再定。
检查点输出也用C语言的二进制IO即可。我一般每2000步把子域数组写成一个独立文件，文件名带rank，用fopen("ckpt_%d.bin", "wb")，配合fwrite一次性写一块连续数组。这里要特别注意：写ghost区域时不要多写，否则多个进程会重复覆盖全局边界数据。读取检查点时进程数必须和写入时相同，或者写一个带头信息的小文件记录全局尺寸和rank数。用文本printf写几十万格点会非常慢，时间步进中尽量只做fwrite。
4. 并行FDTD参数设置与性能优化：从CFL到缓存
4.1 网格步长、CFL和PML厚度之间怎么配合
并行FDTD参数分成物理参数和并行参数。物理参数先按目标频率定网格步长，再按CFL条件定时间步长，最后定PML。如果是在并行环境下调试不稳定，先别急着查代码，先检查是否用了过大的CFLN。PML层数在并行时影响更隐蔽：PML只有在全局外边界才有，内部子域不应再分配PML厚度。一个可取的实现是每个子域在逻辑上都留ghost层，但PML系数数组只在coords等于0或dims-1的进程上分配；否则10层PML会让内存多出约20%到30%，对通信量也是额外负担。
4.2 进程数与本地网格规模：先估内存，再跑起来
确定进程数前先估算每个进程的内存。一个三维并行FDTD程序，每格点需要6个双精度场分量，加上PML辅助数组、材料数组，按每个场分量一套double算，单个子域的场内存在6*8*N_local字节；再加上材料参数和检查点缓冲区，通常按每格点150到200字节估。因此一个进程的N_local建议控制在几十万到几百万，太小通信占比高，太大单步延迟高。编译和运行命令通常长这样：
mpicc -O2 -std=c11 -fopenmp -march=native -o fdtd_mpi fdtd_mpi.c -lm mpirun -np 8 --bind-to core ./fdtd_mpi --grid 4096 2048 --steps 1000-march=native让编译器使用当前CPU的自动向量化扩展，FDTD更新循环对这种优化非常敏感。--bind-to core把每个MPI进程绑到固定核，避免运行中迁移导致缓存击穿。--grid后跟的是全局网格尺寸，MPI_Dims_create会把它拆到8个进程上。如果某个子域的格点数比相邻进程多出10%，考虑换一组进程数因子，例如从2×4改成4×2再测一次。
4.3 restrict指针、内存连续与OpenMP混合并行
纯MPI在每个节点内也走通信协议栈，尽管MPI实现会优化共享内存，但多核节点上OpenMP混合并行更划算。混合模式里，每个MPI进程只负责一个子域，但用OpenMP线程共享该子域数组。下面的更新循环加入OpenMP：
#pragma omp parallel for collapse(2) for (int i = 1; i < nx_local - 1; i++) { for (int j = 1; j < ny_local - 1; j++) { int idx = i * ny_local + j; ez[idx] += coef_ez * (hy[idx] - hy[(i - 1) * ny_local + j]); } }collapse(2)把外层i和j合并成一个大循环，编译器更容易生成静态负载均衡。这里i从1开始，hy[(i - 1)*ny_local + j]访问的是当前子域内部或ghost行，避免数组头越界；边界格点由边界更新单独处理。C语言指针运算不会自动检查边界，这是并行FDTD最常见的崩溃来源。建议在调试期打开-fsanitize=address跑一次小规模网格，能快速抓到越界写。关于C语言指针，我有一条习惯：所有数组基地址只在一处计算，子域内部使用相对索引，不在循环里反复做base + i*ld。这样既防止指针别名，也让MPI通信的发送边界地址容易计算。OpenMP线程数由OMP_NUM_THREADS控制，运行前设OMP_PROC_BIND=true，否则线程迁移可能把缓存局部性毁掉。
4.4 通信与计算重叠的适用场景
3.3节里的Waitall放在更新之前，通信完全串行。对于强可扩展性逼近极限的时候，可以改成：先Isend/Irecv，然后更新子域内部不依赖ghost的格点，再Waitall，最后更新依赖ghost的边界格点。内部区域的循环范围是i从2到nx_local-2，j从2到ny_local-2，避开最外两层；边界格点单独再循环一遍。这个重叠对MPI+OpenMP混合模式尤其有效，因为线程可以在等待通信时继续算内部点。代价是循环被拆成两遍，边界点的重复计算让总浮点次数略增。实测中只有当子域足够大，比如每边超过128格时，重叠的收益才明显；小网格反而因为循环分块损失性能。进入这个阶段前，先用perf简单看下通信时间占比，再决定要不要做。
5. 并行FDTD的验证技巧：一致性、弱可扩展性与调试
5.1 先验证并行一致性，再谈物理正确性
并行FDTD最容易出的是边界索引错乱，而不是物理模型错。所以第一个验证是并行一致性：同一全局网格，分别用1个和4个进程跑相同步数，取某个探针点比较。1个进程的探针值序列是参考，4个进程的任何偏差都说明通信或守护区有问题。实现时可以在每个进程里保存探针坐标落在本子域时的值，最后用MPI_Gather汇总到rank 0输出。如果偏差从某一时间步开始，错的是数据交换；如果一上来就不同，错的是初始场或者探针坐标映射。
5.2 和解析解比较，确认二阶精度没被破坏
并行一致性通过后再和解析解比较，比如平面波穿过均匀介质的衰减，或者点源在远场的1/r衰减。误差一般定义为|E_sim - E_ana|/|E_ana|，在PML边界附近误差可能被边界反射污染，所以探针放在计算域中心区域。如果误差随网格步长减半而接近原来的1/4，说明二阶精度没有因为并行分解破坏；误差不变则多半是PML参数没调好或激励源加载位置不对。
5.3 弱可扩展性测试：保持每进程网格数不变
保持每个MPI进程的子域网格固定，让进程数翻倍，全局网格也随之翻倍，测量单步时间。理想情况下墙钟时间保持不变。下面是一个简易脚本：
for np in 1 2 4 8; do mpirun -np $np --bind-to core ./fdtd_mpi \ --local-grid 128 128 --steps 500 done--local-grid是进程内的局部网格，总网格格点数为128*128*np，单维长度近似128*sqrt(np)。如果平均每步时间随np明显上升，检查是不是某些子域因为全局边界带PML而导致负载不均，或者MPI_Cart的进程排列让跨节点通信路径变长。弱可扩展性的目标不是追求加速比线性，而是单步时间增长控制在20%以内。调试阶段如果出现随进程数增多而偶发的结果漂移，优先怀疑非阻塞通信的发送缓冲区：MPI_Isend在MPI_Waitall之前不能改写send_buf，一旦下个时间步的循环提前覆盖了它，就会出现只有多进程才复现的随机错误。
本文还有配套的精品资源，点击获取
