简介面向计算电磁学与并行编程学习者的 FDTD 方法 C 语言实现以有限差分时域法为核心对麦克斯韦方程组做时间和空间离散适用于天线辐射、点源激励等典型电磁仿真场景尤其适合希望在大规模网格仿真中借助多核并行提升效率的开发者与研究者参考。压缩包总共包含 26 个文件以 h 头文件、cpp 源码、txt 模型与配置文档为主体辅以 docx 使用说明、makefile 构建文件、bat 清理脚本以及可执行程序 MedFDTD.exe包体大小仅 246KB结构清晰便于按需学习。该资源目前已吸引 439 人查看学习。其中内置了点源设置示例用户可在 model.txt 中调整位置、频率与极化参数结合源代码阅读能够完整走通 FDTD 时间步进更新、边界条件处理与并行化改造流程对于想深入研究计算电磁学并行算法或快速搭建电磁仿真验证环境的读者是一份实用且轻量的参考资料。1. 有限差分时域法并行代码先想清楚要不要上MPI再动键盘300×300×300 的网格、2000 个时间步单核 C 语言写的有限差分时域法FDTD程序要跑到第二天换成 8 进程的 MPI 并行版一个上午就能出结果。加速的底气来自数据依赖每个网格点的电场、磁场更新只碰相邻几个点把计算域沿空间轴切开分给不同进程边界上交换一层壳数据即可。这篇讲计算电磁学里最常用的工程路线C 语言 MPI Yee 网格覆盖空间域分解、halo 交换、非阻塞通信和正确性验证。适合已经跑通单线程 FDTD、准备上多进程的工程师也适合正要评估「要不要为这个算例重写并行层」的人。动手前先确认一件事网格总内存超过单机物理内存的四分之一或者单步计算时间远大于一次边界通信的时间再上并行否则先优化缓存局部性更划算。2. 有限差分时域法的Yee网格与空间域分解数据依赖决定切法2.1 Yee网格为什么把电场和磁场错开半个格点有限差分时域法把 Maxwell 旋度方程在时间和空间上做中心差分得到显式递推。1966 年 Yee 提出的网格把 E 和 H 在空间中错开半个格点、在时间上错开半个步长。以真空中的 Ez 分量为例更新式是Ez^{n1}(i,j,k) Ez^n(i,j,k) (dt/eps0)·[ (Hy^n(i1/2,j,k) − Hy^n(i−1/2,j,k))/dx − (Hx^n(i,j1/2,k) − Hx^n(i,j−1/2,k))/dy ]右侧只出现第 n 步的场值下标只偏离半个格点。六个场分量共享同一规律更新一个点只需要半径为一个格点的邻域值。这带来两个直接结论单步计算是纯最近邻操作缓存命中率高单核性能容易优化计算域可以沿任意空间轴切开只有切口两侧的进程需要通信。时间推进是全局同步的dt 受 CFL 条件约束所有进程必须用同一个 dt 推进相同步数所以并行策略只能落在空间切分上时间轴没有可拆的余地。2.2 slab、pencil、cube 三种切分的通信量对比空间域分解的三种基本形态沿一个轴切板slab、沿两个轴切条pencil、沿三个轴切块cube。切分维数越高单块外表面越小总通信量越低代价是邻居关系和索引换算变复杂。总网格 Nx×Ny×Nz、进程数 P 时slab 切分每步每场分量要交换约 2×(Ny×Nz) 个 halo 格点cube 切分时每个方向的面都要交换总量约 2×(Ny×Nz/P^{2/3} Nx×Nz/P^{2/3} Nx×Ny/P^{2/3})。切分方式邻居数512³ 网格、64 进程时每步每场通信量适用规模slab一维2约 4.2M 格点18 进程pencil二维4约 2.2M 格点864 进程cube三维6约 0.1M 格点64 进程以上注意通信量随切分数下降很快但邻居数从 2 涨到 6MPI 调用条数也成倍增加。经验规则是进程数每翻一倍优先考虑增加切分维数而不是继续压缩单维块长后者会同时抬高通信占比和消息条数。2.3 负载均衡按单元数切而不是按坐标均分均匀网格下 NX/nprocs 整数除法够用但天线馈电点附近、介质填充区域的计算量并不均匀。常见做法是先把网格按切分方向做单元数前缀和再按累计单元数找切分点/* 每个格点的计算权重 w[i]按累计单元数找切分点 */ long acc 0, target total_cells / nprocs; int cut 0; for (int i 0; i NX; i) { acc w[i]; if (acc (rank 1) * target) { cut i; break; } }逻辑说明target 是每个进程应分到的累计权重rank 越靠后切点越靠右循环找到第一个超过阈值的格点作为本进程右边界。切分点变化会让各块长度不同halo 交换时收发长度必须按邻居真实尺寸来定所以建议把「切分算法」和「通信代码」解耦切分只输出局部起点和长度通信代码只认局部下标。切完用MPI_Allreduce统计各进程局部单元数最大值与最小值比值超过 1.1 就认为负载不均衡优先调整切分点而不是改通信代码。3. 用C语言写FDTD的MPI并行骨架从域分配到halo交换3.1 用 MPI_Cart_create 把进程组织成三维网格自己算邻居下标容易错常见做法是创建笛卡尔虚拟拓扑让 MPI 接管邻居管理#include stdio.h #include stdlib.h #include mpi.h #define NX 512 #define NY 256 #define NZ 256 int main(int argc, char **argv) { int rank, nprocs; MPI_Init(argc, argv); MPI_Comm_rank(MPI_COMM_WORLD, rank); MPI_Comm_size(MPI_COMM_WORLD, nprocs); int dims[3] {1, 1, 0}; /* 0 表示让库自动分配 */ MPI_Dims_create(nprocs, 3, dims); int periods[3] {0, 0, 0}; /* 散射问题全部不周期 */ MPI_Comm cart; MPI_Cart_create(MPI_COMM_WORLD, 3, dims, periods, 0, cart); int coords[3]; MPI_Cart_coords(cart, rank, 3, coords); printf(rank %d - coords (%d,%d,%d)\n, rank, coords[0], coords[1], coords[2]); /* 后续时间步推进全部用 cart 通信子 */ ... }参数说明dims 数组里写 0 的位置由MPI_Dims_create自动补它会尽量让三个维度均衡避免出现某维只有 1 个进程的畸形拓扑periods 各方向是否开周期要和物理边界一致波导模拟沿传播方向置 1散射计算四周留给 CPML 必须全 0。MPI_Cart_create之后不要再拿MPI_COMM_WORLD下发通信调用否则消息会同其他组串扰。3.2 局部数组与全局下标换算沿 x 方向一维切分每个进程持有 (nx2)×(ny2)×(nz2) 的数组含两层 halo/* 沿 x 方向切分每个进程持有 (local_nx2)*(NY2)*(NZ2) 数组 */ int local_nx NX / nprocs; int rem NX % nprocs; if (rank rem) local_nx; /* 前 rem 个进程多拿一层 */ int global_start 0; for (int p 0; p rank; p) global_start NX / nprocs (p rem ? 1 : 0); double *Ez malloc((local_nx 2) * (NY 2) * (NZ 2) * sizeof(double)); /* 全局下标 i_global global_start (i_local - 1)减 1 是因为 i_local0 是 halo */说明局部下标 0 和 local_nx1 是 halo内部格点是 1..local_nx所以全局换算要减 1。NX 不能被进程数整除时各块长度不同halo 交换的发送长度必须按邻居实际尺寸计算。初始化时把各进程的 local_nx 用MPI_Allgather收集一份构造收发缓冲时直接用邻居值不要从自己的 local_nx 去推。这里的 malloc 就是 C 语言内存管理的重点局部数组反复 malloc/free 会带来页错误抖动时间步循环外只分配一次循环内复用。3.3 halo 交换的 MPI_Sendrecv 实现Yee 网格里 x 方向的邻居交换关键是 yz 面在内存中不连续要用派生类型#define IDX(i,j,k) ((i)*(NY2)*(NZ2) (j)*(NZ2) (k)) /* 交换 Ez 在 x 方向的左右 halo */ void exchange_halo_x(double *Ez, int nx, int ny, int nz, MPI_Comm cart) { int left, right; MPI_Cart_shift(cart, 0, 1, left, right); static MPI_Datatype face MPI_DATATYPE_NULL; if (face MPI_DATATYPE_NULL) { /* ny 行每行 nz 个双精度行间距 nz2 */ MPI_Type_vector(ny, nz, nz 2, MPI_DOUBLE, face); MPI_Type_commit(face); } int tag 10; /* 把 i1 内部面发给 left从 left 收到数据放进 i0 的 halo */ MPI_Sendrecv(Ez[IDX(1,1,1)], 1, face, left, tag, Ez[IDX(0,1,1)], 1, face, left, tag, cart, MPI_STATUS_IGNORE); /* 把 inx 内部面发给 right从 right 收到数据放进 inx1 的 halo */ MPI_Sendrecv(Ez[IDX(nx,1,1)], 1, face, right, tag, Ez[IDX(nx1,1,1)], 1, face, right, tag, cart, MPI_STATUS_IGNORE); }代码说明MPI_Sendrecv成对收发避免先 Send 后 Recv 在大消息下死锁。MPI_Type_vector用「指针 步长」描述不连续面一次调用收发整个面省掉临时缓冲的整层拷贝。tag 必须按方向区分否则相邻两个方向的同字段消息会错位工程上我给 x/y/z 三个方向分别用 tag10/11/12字段名写进调试开关。3.4 主时间步循环里「先场更新、后 halo 交换」for (int t 0; t nsteps; t) { /* 第 n1/2 步H 由上一轮已就绪的 E 更新只用本进程内部点 */ update_H(Hx, Hy, Hz, Ex, Ey, Ez, dt, dx, dy, dz); exchange_halo_x(Hx); exchange_halo_y(Hx); exchange_halo_z(Hx); /* Hy、Hz 同样处理 */ /* 第 n1 步E 由刚交换完的 H 更新 */ update_E(Ex, Ey, Ez, Hx, Hy, Hz, dt, dx, dy, dz); exchange_halo_x(Ex); exchange_halo_y(Ex); exchange_halo_z(Ex); /* Ey、Ez 同样处理 */ if (t % 100 0 rank 0) printf(step %d done\n, t); }顺序不能调换H 更新依赖的是上一轮结束前已经交换好的 E 的 halo如果先更新 H 再交换 E切口处会差出整整一轮的数据。每个场分量算完立刻交换不要等六个场分量全部更新完再统一通信那会让在途消息数翻倍64 进程以上时明显损伤链路利用率。dt 由全局最小网格间距决定初始化时用MPI_Allreduce求全局最小 dx、dy、dz 后统一广播避免各进程自行计算出现不一致。4. 并行FDTD的关键参数通信缓冲、非阻塞收发与CPML边界4.1 通信量随切分维数上升而下降的定量规律每步每场分量的通信量近似等于 2 × halo 层数 × 单块切面面积。切成 (px, py, pz) 块时通信面正比于 (Ny×Nz)/px (Nx×Nz)/py (Nx×Ny)/pz。直观理解切得越碎单块外表面越小但切面总数变多总通信面积变化不大真正决定通信时间的是「单条消息的大小 × 消息条数」。slab 切分单条消息大但条数少cube 切分条数多但每条小MPI 的短消息延迟在百微秒量级所以进程数上去后必须降低消息条数。判据很简单当单步计算时间降到与两倍消息延迟同量级时就该增加切分维数。4.2 用 MPI_Isend/Irecv 把通信藏进场更新halo 交换的数据只影响下几步的边界格点内部格点的更新完全不依赖邻居这给计算通信重叠留了空间MPI_Request reqs[2]; int tag 20; /* 先发起异步收发不等待 */ MPI_Isend(Ez[IDX(nx,1,1)], 1, face, right, tag, cart, reqs[0]); MPI_Irecv(Ez[IDX(nx1,1,1)], 1, face, right, tag, cart, reqs[1]); /* 中间插入与 x 方向 halo 无关的计算y/z 方向的内部点更新 */ update_E_inner_yz(Ex, Ey, Ez, Hx, Hy, Hz, dt); MPI_Waitall(2, reqs, MPI_STATUSES_IGNORE);说明Isend/Irecv 配对后MPI_Waitall成对等待。注意一个经典坑MPI_Isend之后如果不 Wait 就释放发送缓冲区小消息常被 MPI 内部缓冲而「成功」大消息可能直接走同步协议导致死锁。缓冲区生命周期必须覆盖到 Waitall 之后。重叠收益在 slab 切分、单条消息大时最明显cube 切分后消息变小重叠收益下降这时代价是代码里多了 reqs 数组管理建议先用阻塞版跑通正确性再改非阻塞版。4.3 CPML 吸收边界在并行域里的归属并行 FDTD 必须有吸收边界常见做法是 CPML。CPML 每个方向要额外维护 psi 辅助数组内存开销比真空区域高约一倍。域分解时 CPML 区域要当作普通网格参与切分不能单独开进程——否则 CPML 那几个进程负载显著低于内部区域破坏负载均衡。CPML 辅助场的更新不涉及跨进程项只有主 E/H 场需要 halo所以它只影响内存预算不影响通信代码结构。网格四个侧面都是 CPML 时角点格点同时属于两个方向的 CPML辅助数组要各算各的更新公式里叠加两个方向的 psi 贡献。4.4 并行FDTD的参数表参数推荐取值说明halo 层数1 层调试期可设 2 层1 层满足 Yee 最近邻依赖2 层便于一致性校验每步通信次数6 个场分量各 1 次超过 6 次说明消息拆碎了MPI_DatatypeType_vector 构造面类型避免整层拷入临时缓冲切分维数≤8 进程一维864 二维64 三维参考 2.2 通信量公式时间步内 Barrier仅调试用进程多时 Barrier 本身成为瓶颈消息 tag按方向分配 10/11/12防止相邻方向消息错位5. 用探针点电场和全局能量守恒验证并行FDTD没写错5.1 单进程参照解是并行代码的标尺并行版本最容易错在下标换算和 halo 收发方向。先用单进程跑一个点源算例记录某个探针点的 Ez 时间序列作为基准再用np2、np4跑相同算例探针点如果落在某个进程内部直接对比该点输出。差异超过 1e-10 说明切分或通信有错不需要看完整场分布就能定位。5.2 用全局能量守恒在每 100 步做一次自检无耗散真空里总能量应当恒定这是校验 halo 是否丢数据的最灵敏指标double E_energy 0.0, H_energy 0.0; for (int i 1; i nx; i) for (int j 1; j ny; j) for (int k 1; k nz; k) { E_energy eps0 * (Ex*Ex Ey*Ey Ez*Ez); H_energy mu0 * (Hx*Hx Hy*Hy Hz*Hz); } double local 0.5 * (E_energy H_energy); /* E、H 半格时间交错 */ double total; MPI_Allreduce(local, total, 1, MPI_DOUBLE, MPI_SUM, cart); if (rank 0 t % 100 0) printf(step %d energy %.12e\n, t, total);参数说明MPI_Allreduce用 MPI_SUM 累加各进程能量任何 halo 丢失都会让总能量出现阶跃式下降。真空里能量应恒定到小数点后 8 位以上加了 CPML 后能量单调下降突然跳变说明通信错位。5.3 加速比测试与 halo 一致性断言跑强扩展测试时固定网格规模用time mpirun -np 4 ./fdtd这类命令记录每千步墙钟时间完整算例至少跑 1000 步再计时避免进程启动抖动。加速比接近进程数说明通信占比低加速比掉到 0.7 以下先看消息条数再看是否还有 Barrier 残留。最后一个技巧在调试版本里加一段 halo 一致性断言——进程 A 发出去的内部面进程 B 收到后放进 halo下一轮 A 再把自己的内部面重发一次用MPI_Compare_and_swap或直接在 B 侧逐元素对比两个来源的差异超过 1e-12 就MPI_Abort。这个断言能瞬间区分「切分点算错」和「收发方向接反」两类最隐蔽的错误等完全跑稳再关掉编译开关。本文还有配套的精品资源点击获取