高频电磁场仿真-主题072-计算电磁学高性能计算

引言

　　计算电磁学（Computational Electromagnetics, CEM）是电磁场理论、数值方法与计算机科学交叉的前沿领域。随着电磁系统复杂度的增加和精度要求的提高，传统串行计算已难以满足实际工程需求。

计算电磁学的计算挑战

　　典型计算规模：

• 电小尺寸问题：波长λ = 1m，目标尺寸0.1λ，网格尺寸λ/20，总网格数约800个
• 电中等问题：目标尺寸5λ，网格数约10^6个
• 电大问题：飞机雷达散射（10m长，3GHz），网格数约10^9个
• 超电大问题：城市环境传播，网格数可达10^12个

　　计算复杂度分析：

• FDTD方法：O(N)每时间步，总复杂度O(N × N_t)
• 矩量法：O(N²)内存，O(N³)或O(N² log N)求解
• 有限元法：O(N)内存，O(N^1.5)求解（直接法）

高性能计算发展历程

1970s-1980s：向量计算机（Cray-1）
    ↓
1990s：MPP大规模并行机（IBM SP2）
    ↓
2000s：集群计算（Beowulf集群）
    ↓
2010s：多核CPU + GPU异构计算
    ↓
2020s：E级计算（Exascale）、AI加速器

性能指标与评估

　　关键性能指标：

• 峰值性能（Peak Performance）：理论最大计算能力（FLOPS）
• 持续性能（Sustained Performance）：实际应用中的性能
• 并行效率（Parallel Efficiency）：
  η = T₁ / (p × Tₚ)
• 加速比（Speedup）：
  S(p) = T₁ / Tₚ

　　Amdahl定律：

S(p) = 1 / ((1−f) + f/p)

其中f为可并行化比例，p为处理器数量。

高性能计算基础

计算机体系结构

　　现代处理器架构：

1. 多核CPU架构：
   • 共享缓存层次（L1/L2/L3）
   • SIMD向量单元（AVX-512, NEON）
   • 超标量流水线与乱序执行
2. 内存层次结构：
   

   寄存器（1 cycle）
       ↓
   L1缓存（4 cycles）
       ↓
   L2缓存（10 cycles）
       ↓
   L3缓存（40 cycles）
       ↓
   主内存（200 cycles）
       ↓
   磁盘（10^6 cycles）
   

3. GPU架构：
   • 大量轻量级线程（数千个）
   • 高带宽显存（HBM2e/HBM3）
   • SIMT执行模型

并行计算模型

类型	描述	示例
SISD	单指令单数据	传统串行处理器
SIMD	单指令多数据	GPU、向量处理器
MISD	多指令单数据	容错系统（少见）
MIMD	多指令多数据	多核CPU、集群

　　并行编程模型：

1. 共享内存模型（OpenMP）：
   • 线程级并行
   • 统一地址空间
   • 适合多核CPU
2. 分布式内存模型（MPI）：
   • 进程级并行
   • 显式消息传递
   • 适合集群系统
3. 异构编程模型（CUDA/OpenCL）：
   • 主机+设备架构
   • 显存管理
   • 适合GPU加速

性能优化原则

　　Roofline模型：

Performance = min{Peak FLOPS, Bandwidth × Operational Intensity}

　　优化策略层次：

1. 算法级优化：
   • 选择合适算法（O(N) vs O(N³)）
   • 快速算法（FFT、FMM）
2. 实现级优化：
   • 循环优化（向量化、展开）
   • 内存访问优化（局部性、对齐）
3. 并行级优化：
   • 负载均衡
   • 通信优化
   • 同步最小化

并行计算模型与架构

域分解方法

　　空间域分解（Domain Decomposition）：

对于FDTD等基于网格的方法，将计算域划分为多个子域：

┌─────────┬─────────┬─────────┐
│  Sub0   │  Sub1   │  Sub2   │
│ (0:Nx/3)│(Nx/3:2Nx/3)│(2Nx/3:Nx)│
├─────────┼─────────┼─────────┤
│  Sub3   │  Sub4   │  Sub5   │
├─────────┼─────────┼─────────┤
│  Sub6   │  Sub7   │  Sub8   │
└─────────┴─────────┴─────────┘

　　分解策略：

• 一维分解：沿x轴分割，适合细长线状结构
• 二维分解：xy平面分割，适合面状结构
• 三维分解：xyz立体分割，适合体状结构

　　边界数据交换：

每个子域需要与相邻子域交换边界层数据（halo exchange）：

# 伪代码：边界交换
for step in range(n_steps):
    # 计算内部点
    compute_interior(subdomain)
    
    # 非阻塞发送边界数据
    send_halo_to_neighbors(subdomain)
    
    # 接收边界数据
    recv_halo_from_neighbors(subdomain)
    
    # 计算边界点
    compute_boundary(subdomain)

负载均衡

　　静态负载均衡：

• 均匀网格划分：每个进程相同网格数
• 非均匀划分：根据计算复杂度调整

　　动态负载均衡：

• 自适应网格细化（AMR）
• 工作窃取（Work Stealing）

　　负载均衡度量：

L = maxᵢ(Tᵢ) / T̄

其中Tᵢ为进程i的计算时间，理想情况下L ≈ 1。

通信优化

　　通信模式：

1. 点对点通信：
   • MPI_Send / MPI_Recv
   • 阻塞与非阻塞
2. 集合通信：
   • MPI_Bcast：广播
   • MPI_Reduce：归约
   • MPI_Allgather：全收集
   • MPI_Alltoall：全交换

　　通信优化技术：

1. 通信隐藏：
   

   # 计算与通信重叠
   compute_interior()      # 计算内部（不依赖边界）
   start_halo_exchange()   # 启动非阻塞通信
   compute_boundary()      # 计算边界（需要等待通信完成）
   wait_halo_exchange()    # 等待通信完成
   

2. 聚合通信：
   • 减少小消息数量
   • 使用打包/解包
3. 拓扑感知映射：
   • 将相邻子域映射到相邻节点
   • 减少网络跳数

OpenMP并行FDTD

OpenMP基础

　　OpenMP编程模型：

OpenMP是一种基于编译指令的共享内存并行编程接口。

#pragma omp parallel for
for (int i = 0; i < N; i++) {
    // 并行执行的循环体
    result[i] = compute(data[i]);
}

　　关键指令：

• #pragma omp parallel：创建并行区域
• #pragma omp for：循环并行化
• #pragma omp sections：任务并行
• #pragma omp critical：临界区
• #pragma omp barrier：同步点

一维FDTD并行化

　　算法原理：

import numpy as np
from numba import njit, prange
import time

@njit(parallel=True)
def fdtd_1d_parallel(E, H, ce, ch, n_steps, nx):
    """
    一维FDTD并行实现（使用Numba的prange）
    """
    for n in range(n_steps):
        # 并行更新磁场
        for i in prange(nx - 1):
            H[i] = H[i] + ch[i] * (E[i + 1] - E[i])
        
        # 并行更新电场
        for i in prange(1, nx):
            E[i] = E[i] + ce[i] * (H[i] - H[i - 1])
        
        # 激励源（在中心位置）
        E[nx // 2] += np.sin(2 * np.pi * 1e9 * n * 1e-12)
    
    return E, H

# 参数设置
nx = 10000
n_steps = 1000
E = np.zeros(nx)
H = np.zeros(nx - 1)
ce = np.ones(nx) * 3e8 * 1e-12 / 1e-3  # 光速×时间步/空间步
ch = np.ones(nx - 1) * 3e8 * 1e-12 / 1e-3

# 运行并行FDTD
start = time.time()
E_result, H_result = fdtd_1d_parallel(E, H, ce, ch, n_steps, nx)
parallel_time = time.time() - start

print(f"并行计算时间: {parallel_time:.3f}秒")

三维FDTD并行化

　　三维网格更新：

@njit(parallel=True)
def update_H_3d_parallel(Hx, Hy, Hz, Ex, Ey, Ez, ch, nx, ny, nz):
    """
    并行更新三维磁场分量
    """
    # Hx更新
    for i in prange(nx):
        for j in range(ny - 1):
            for k in range(nz - 1):
                Hx[i, j, k] += ch * ((Ez[i, j + 1, k] - Ez[i, j, k]) - 
                                     (Ey[i, j, k + 1] - Ey[i, j, k]))
    
    # Hy更新
    for i in prange(nx - 1):
        for j in range(ny):
            for k in range(nz - 1):
                Hy[i, j, k] += ch * ((Ex[i, j, k + 1] - Ex[i, j, k]) - 
                                     (Ez[i + 1, j, k] - Ez[i, j, k]))
    
    # Hz更新
    for i in prange(nx - 1):
        for j in range(ny - 1):
            for k in range(nz):
                Hz[i, j, k] += ch * ((Ey[i + 1, j, k] - Ey[i, j, k]) - 
                                     (Ex[i, j + 1, k] - Ex[i, j, k]))

@njit(parallel=True)
def update_E_3d_parallel(Ex, Ey, Ez, Hx, Hy, Hz, ce, nx, ny, nz):
    """
    并行更新三维电场分量
    """
    # Ex更新
    for i in prange(nx):
        for j in range(1, ny):
            for k in range(1, nz):
                Ex[i, j, k] += ce * ((Hz[i, j, k] - Hz[i, j - 1, k]) - 
                                     (Hy[i, j, k] - Hy[i, j, k - 1]))
    
    # Ey更新
    for i in prange(1, nx):
        for j in range(ny):
            for k in range(1, nz):
                Ey[i, j, k] += ce * ((Hx[i, j, k] - Hx[i, j, k - 1]) - 
                                     (Hz[i, j, k] - Hz[i - 1, j, k]))
    
    # Ez更新
    for i in prange(1, nx):
        for j in range(1, ny):
            for k in range(nz):
                Ez[i, j, k] += ce * ((Hy[i, j, k] - Hy[i - 1, j, k]) - 
                                     (Hx[i, j, k] - Hx[i, j - 1, k]))

向量化优化

　　SIMD向量化：

现代CPU支持AVX-512（512位向量），可同时处理8个双精度或16个单精度浮点数。

import numpy as np
from numba import njit, vectorize, float64

@vectorize([float64(float64, float64, float64)], target='parallel')
def vectorized_update(E, H_diff, coeff):
    """
    向量化电场更新
    """
    return E + coeff * H_diff

# 使用示例
E_new = vectorized_update(E, H[1:] - H[:-1], ce)

MPI分布式计算

MPI基础

　　MPI程序结构：

from mpi4py import MPI

comm = MPI.COMM_WORLD
rank = comm.Get_rank()
size = comm.Get_size()

print(f"Process {rank} of {size}")

# 并行计算...

comm.Barrier()  # 同步
MPI.Finalize()

　　基本通信操作：

# 点对点通信
if rank == 0:
    data = np.array([1.0, 2.0, 3.0])
    comm.Send(data, dest=1, tag=0)
elif rank == 1:
    data = np.empty(3, dtype=np.float64)
    comm.Recv(data, source=0, tag=0)

# 非阻塞通信
req = comm.Isend(data, dest=1, tag=0)
# ... 做其他计算 ...
req.Wait()

# 广播
if rank == 0:
    data = np.array([1.0, 2.0, 3.0])
else:
    data = np.empty(3, dtype=np.float64)
comm.Bcast(data, root=0)

# 归约
local_sum = np.array([rank], dtype=np.float64)
global_sum = np.empty(1, dtype=np.float64)
comm.Allreduce(local_sum, global_sum, op=MPI.SUM)

一维FDTD MPI实现

　　域分解与通信：

from mpi4py import MPI
import numpy as np

class FDTD1D_MPI:
    def __init__(self, nx_global, n_steps):
        self.comm = MPI.COMM_WORLD
        self.rank = self.comm.Get_rank()
        self.size = self.comm.Get_size()
        
        self.nx_global = nx_global
        self.n_steps = n_steps
        
        # 计算本地网格范围
        self.nx_local = nx_global // self.size
        self.start_idx = self.rank * self.nx_local
        self.end_idx = self.start_idx + self.nx_local
        
        # 添加幽灵单元（halo）
        self.E = np.zeros(self.nx_local + 2)  # +2用于边界交换
        self.H = np.zeros(self.nx_local + 2)
        
        self.ce = 0.5
        self.ch = 0.5
    
    def exchange_boundaries(self):
        """
        与相邻进程交换边界数据
        """
        # 发送左边界到左邻居，接收右边界来自右邻居
        if self.rank > 0:
            self.comm.Send(self.E[1:2], dest=self.rank-1, tag=0)
            self.comm.Recv(self.E[0:1], source=self.rank-1, tag=1)
        
        # 发送右边界到右邻居，接收左边界来自左邻居
        if self.rank < self.size - 1:
            self.comm.Send(self.E[-2:-1], dest=self.rank+1, tag=1)
            self.comm.Recv(self.E[-1:], source=self.rank+1, tag=0)
    
    def update_fields(self):
        """
        更新电磁场
        """
        # 更新H（内部点）
        for i in range(1, self.nx_local + 1):
            self.H[i] += self.ch * (self.E[i+1] - self.E[i])
        
        # 更新E（内部点）
        for i in range(1, self.nx_local + 1):
            self.E[i] += self.ce * (self.H[i] - self.H[i-1])
    
    def add_source(self, n):
        """
        在中央位置添加源
        """
        center = self.nx_global // 2
        if self.start_idx <= center < self.end_idx:
            local_idx = center - self.start_idx + 1
            self.E[local_idx] += np.sin(2 * np.pi * n * 0.1)
    
    def run(self):
        """
        运行FDTD仿真
        """
        for n in range(self.n_steps):
            self.exchange_boundaries()
            self.update_fields()
            self.add_source(n)
        
        return self.E[1:-1]  # 返回内部点

# 主程序
if __name__ == '__main__':
    fdtd = FDTD1D_MPI(nx_global=10000, n_steps=1000)
    result = fdtd.run()
    
    if fdtd.rank == 0:
        print(f"MPI FDTD completed on {fdtd.size} processes")

三维FDTD MPI实现

　　Cartesian拓扑：

from mpi4py import MPI
import numpy as np

class FDTD3D_MPI:
    def __init__(self, nx, ny, nz, n_steps, px=2, py=2, pz=2):
        """
        3D FDTD with MPI Cartesian topology
        
        px, py, pz: 各维度的进程数
        """
        self.comm = MPI.COMM_WORLD
        self.rank = self.comm.Get_rank()
        self.size = self.comm.Get_size()
        
        assert self.size == px * py * pz
        
        self.nx, self.ny, self.nz = nx, ny, nz
        self.n_steps = n_steps
        
        # 创建3D笛卡尔拓扑
        self.dims = [px, py, pz]
        self.periods = [False, False, False]
        self.cart_comm = self.comm.Create_cart(self.dims, self.periods, reorder=True)
        
        self.coords = self.cart_comm.Get_coords(self.rank)
        
        # 计算本地网格大小
        self.nx_local = nx // px
        self.ny_local = ny // py
        self.nz_local = nz // pz
        
        # 分配场数组（含halo）
        self.Ex = np.zeros((self.nx_local + 2, self.ny_local + 2, self.nz_local + 2))
        self.Ey = np.zeros((self.nx_local + 2, self.ny_local + 2, self.nz_local + 2))
        self.Ez = np.zeros((self.nx_local + 2, self.ny_local + 2, self.nz_local + 2))
        
        self.Hx = np.zeros((self.nx_local + 2, self.ny_local + 2, self.nz_local + 2))
        self.Hy = np.zeros((self.nx_local + 2, self.ny_local + 2, self.nz_local + 2))
        self.Hz = np.zeros((self.nx_local + 2, self.ny_local + 2, self.nz_local + 2))
    
    def exchange_halo(self):
        """
        6个方向的halo交换
        """
        # x方向
        left, right = self.cart_comm.Shift(0, 1)
        
        # 发送右边界到右邻居，从左邻居接收左边界
        if right != MPI.PROC_NULL:
            self.cart_comm.Send(self.Ex[-2:-1, :, :], dest=right, tag=0)
        if left != MPI.PROC_NULL:
            self.cart_comm.Recv(self.Ex[0:1, :, :], source=left, tag=0)
        
        # y方向
        down, up = self.cart_comm.Shift(1, 1)
        # ... 类似处理 ...
        
        # z方向
        back, front = self.cart_comm.Shift(2, 1)
        # ... 类似处理 ...
    
    def update_H(self):
        """
        更新磁场
        """
        nx, ny, nz = self.nx_local, self.ny_local, self.nz_local
        
        # Hx
        self.Hx[1:nx+1, 0:ny, 0:nz] += 0.5 * (
            (self.Ez[1:nx+1, 1:ny+1, 0:nz] - self.Ez[1:nx+1, 0:ny, 0:nz]) -
            (self.Ey[1:nx+1, 0:ny, 1:nz+1] - self.Ey[1:nx+1, 0:ny, 0:nz])
        )
        
        # Hy, Hz类似...
    
    def update_E(self):
        """
        更新电场
        """
        nx, ny, nz = self.nx_local, self.ny_local, self.nz_local
        
        # Ex
        self.Ex[0:nx, 1:ny+1, 1:nz+1] += 0.5 * (
            (self.Hz[0:nx, 1:ny+1, 1:nz+1] - self.Hz[0:nx, 0:ny, 1:nz+1]) -
            (self.Hy[0:nx, 1:ny+1, 1:nz+1] - self.Hy[0:nx, 1:ny+1, 0:nz])
        )
        
        # Ey, Ez类似...
    
    def run(self):
        """
        主循环
        """
        for n in range(self.n_steps):
            self.update_H()
            self.exchange_halo()
            self.update_E()
            
            # 添加源...
            
            if self.rank == 0 and n % 100 == 0:
                print(f"Step {n}/{self.n_steps}")

GPU加速计算

GPU架构与CUDA编程

特性	CPU	GPU
核心数	少（4-64）	多（数千）
线程管理	复杂（乱序执行）	简单（SIMT）
缓存	大（MB级）	小（KB级）
内存带宽	中等（~100 GB/s）	高（~1 TB/s）
适合任务	串行/分支多	并行/计算密集

　　CUDA编程模型：

import cupy as cp

# 创建GPU数组
E_gpu = cp.zeros((1000, 1000), dtype=cp.float32)
H_gpu = cp.zeros((1000, 1000), dtype=cp.float32)

# 定义CUDA核函数
fdtd_kernel = cp.RawKernel(r'''
extern "C" __global__
void fdtd_update(float* E, float* H, float ce, float ch, int nx, int ny) {
    int i = blockDim.x * blockIdx.x + threadIdx.x;
    int j = blockDim.y * blockIdx.y + threadIdx.y;
    
    if (i < nx && j < ny) {
        int idx = i * ny + j;
        
        // 更新H
        if (j < ny - 1) {
            H[idx] += ch * (E[idx + 1] - E[idx]);
        }
        
        // 更新E
        if (i > 0 && j > 0) {
            E[idx] += ce * (H[idx] - H[idx - 1]);
        }
    }
}
''', 'fdtd_update')

# 启动核函数
threads = (16, 16)
blocks = ((nx + 15) // 16, (ny + 15) // 16)
fdtd_kernel(blocks, threads, (E_gpu, H_gpu, ce, ch, nx, ny))

CuPy实现FDTD

　　使用CuPy简化GPU编程：

import cupy as cp
import numpy as np
import time

class FDTD2D_GPU:
    def __init__(self, nx, ny, n_steps):
        self.nx = nx
        self.ny = ny
        self.n_steps = n_steps
        
        # 在GPU上分配数组
        self.Ex = cp.zeros((nx, ny), dtype=cp.float32)
        self.Ey = cp.zeros((nx, ny), dtype=cp.float32)
        self.Hz = cp.zeros((nx, ny), dtype=cp.float32)
        
        self.ce = 0.5
        self.ch = 0.5
    
    def update_H(self):
        """
        更新磁场（使用CuPy的向量化操作）
        """
        self.Hz[:-1, :-1] += self.ch * (
            (self.Ex[:-1, 1:] - self.Ex[:-1, :-1]) -
            (self.Ey[1:, :-1] - self.Ey[:-1, :-1])
        )
    
    def update_E(self):
        """
        更新电场
        """
        self.Ex[1:, :-1] += self.ce * (self.Hz[1:, :-1] - self.Hz[:-1, :-1])
        self.Ey[:-1, 1:] -= self.ce * (self.Hz[:-1, 1:] - self.Hz[:-1, :-1])
    
    def add_source(self, n):
        """
        添加高斯脉冲源
        """
        t0 = 20
        spread = 6
        source = cp.exp(-0.5 * ((n - t0) / spread) ** 2)
        self.Ez[self.nx // 2, self.ny // 2] += source
    
    def run(self):
        """
        主循环
        """
        for n in range(self.n_steps):
            self.update_H()
            self.update_E()
            self.add_source(n)
        
        # 将结果复制回CPU
        return cp.asnumpy(self.Ez)

# 性能对比
def benchmark_gpu():
    sizes = [512, 1024, 2048, 4096]
    
    for nx in sizes:
        print(f"\nGrid size: {nx} x {nx}")
        
        # CPU版本
        fdtd_cpu = FDTD2D_CPU(nx, nx, 100)
        start = time.time()
        result_cpu = fdtd_cpu.run()
        cpu_time = time.time() - start
        print(f"CPU time: {cpu_time:.3f}s")
        
        # GPU版本
        fdtd_gpu = FDTD2D_GPU(nx, nx, 100)
        start = time.time()
        result_gpu = fdtd_gpu.run()
        gpu_time = time.time() - start
        print(f"GPU time: {gpu_time:.3f}s")
        
        speedup = cpu_time / gpu_time
        print(f"Speedup: {speedup:.2f}x")

多GPU并行

　　使用NCCL进行多GPU通信：

import cupy as cp
from cupy.cuda import nccl
import os

class FDTD_MultiGPU:
    def __init__(self, nx, ny, n_steps, n_gpus):
        self.nx = nx
        self.ny = ny
        self.n_steps = n_steps
        self.n_gpus = n_gpus
        
        # 每个GPU处理的部分
        self.ny_local = ny // n_gpus
        
        # 初始化NCCL
        self.nccl_comm = nccl.NcclCommunicator(n_gpus)
        
        self.streams = [cp.cuda.Stream(device=i) for i in range(n_gpus)]
    
    def run(self):
        for n in range(self.n_steps):
            # 每个GPU计算本地部分
            for gpu_id in range(self.n_gpus):
                with cp.cuda.Device(gpu_id):
                    with self.streams[gpu_id]:
                        self.update_fields_gpu(gpu_id)
            
            # 同步
            for stream in self.streams:
                stream.synchronize()
            
            # 交换边界
            self.exchange_boundaries_nccl()

异构计算与负载均衡

CPU+GPU协同计算

　　任务划分策略：

import concurrent.futures
import cupy as cp
import numpy as np

class HeterogeneousFDTD:
    def __init__(self, nx, ny, ratio_gpu=0.7):
        """
        ratio_gpu: GPU处理的网格比例
        """
        self.nx = nx
        self.ny = ny
        self.ratio_gpu = ratio_gpu
        
        # GPU部分
        self.ny_gpu = int(ny * ratio_gpu)
        self.E_gpu = cp.zeros((nx, self.ny_gpu))
        self.H_gpu = cp.zeros((nx, self.ny_gpu))
        
        # CPU部分
        self.ny_cpu = ny - self.ny_gpu
        self.E_cpu = np.zeros((nx, self.ny_cpu))
        self.H_cpu = np.zeros((nx, self.ny_cpu))
    
    def step(self):
        # 启动GPU计算
        gpu_future = self.executor.submit(self.update_gpu)
        
        # CPU并行计算
        self.update_cpu()
        
        # 等待GPU完成
        gpu_future.result()
        
        # 交换边界数据
        self.exchange_cpu_gpu_boundary()
    
    def update_gpu(self):
        with cp.cuda.Device(0):
            # GPU更新代码
            self.H_gpu[:-1, :-1] += 0.5 * (self.E_gpu[:-1, 1:] - self.E_gpu[:-1, :-1])
            # ...
    
    def update_cpu(self):
        # 使用OpenMP并行
        self.H_cpu[:-1, :-1] += 0.5 * (self.E_cpu[:-1, 1:] - self.E_cpu[:-1, :-1])
        # ...

动态负载均衡

　　工作窃取（Work Stealing）：

from queue import Queue
import threading

class WorkStealingScheduler:
    def __init__(self, n_workers):
        self.n_workers = n_workers
        self.local_queues = [Queue() for _ in range(n_workers)]
        self.locks = [threading.Lock() for _ in range(n_workers)]
    
    def submit(self, task, worker_id):
        self.local_queues[worker_id].put(task)
    
    def steal_work(self, thief_id):
        """
        从其他工作线程窃取任务
        """
        for victim_id in range(self.n_workers):
            if victim_id == thief_id:
                continue
            
            with self.locks[victim_id]:
                if not self.local_queues[victim_id].empty():
                    return self.local_queues[victim_id].get()
        
        return None

性能优化技术

内存优化

　　缓存优化策略：

import numpy as np
from numba import njit

@njit(cache=True)
def cache_optimized_update(E, H, nx, ny):
    """
    缓存友好的场更新
    按行访问以提高空间局部性
    """
    for i in range(nx - 1):
        for j in range(ny - 1):
            # 连续内存访问
            H[i, j] += 0.5 * (E[i, j + 1] - E[i, j])

# 分块（Tiling）优化
@njit(parallel=True)
def tiled_update(E, H, nx, ny, tile_size=64):
    """
    分块更新以提高缓存利用率
    """
    for i0 in prange(0, nx, tile_size):
        for j0 in range(0, ny, tile_size):
            # 处理一个块
            i_end = min(i0 + tile_size, nx - 1)
            j_end = min(j0 + tile_size, ny - 1)
            
            for i in range(i0, i_end):
                for j in range(j0, j_end):
                    H[i, j] += 0.5 * (E[i, j + 1] - E[i, j])

通信优化

　　聚合通信与压缩：

import numpy as np
from mpi4py import MPI

class OptimizedCommunication:
    def __init__(self, comm):
        self.comm = comm
        self.rank = comm.Get_rank()
    
    def pack_halo(self, field, halo_width=1):
        """
        打包边界数据
        """
        # 左边界
        left_halo = field[:, :halo_width].copy()
        # 右边界
        right_halo = field[:, -halo_width:].copy()
        
        return left_halo, right_halo
    
    def nonblocking_exchange(self, field):
        """
        非阻塞边界交换
        """
        left_halo, right_halo = self.pack_halo(field)
        
        # 启动非阻塞发送
        requests = []
        
        if self.rank > 0:
            req = self.comm.Isend(left_halo, dest=self.rank-1, tag=0)
            requests.append(req)
        
        if self.rank < self.comm.Get_size() - 1:
            req = self.comm.Isend(right_halo, dest=self.rank+1, tag=1)
            requests.append(req)
        
        # 接收数据
        if self.rank > 0:
            recv_left = np.empty_like(left_halo)
            self.comm.Recv(recv_left, source=self.rank-1, tag=1)
            field[:, 0] = recv_left[:, 0]
        
        if self.rank < self.comm.Get_size() - 1:
            recv_right = np.empty_like(right_halo)
            self.comm.Recv(recv_right, source=self.rank+1, tag=0)
            field[:, -1] = recv_right[:, 0]
        
        # 等待发送完成
        MPI.Request.Waitall(requests)

计算与通信重叠

　　双缓冲技术：

class OverlappingFDTD:
    def __init__(self, nx, ny):
        self.nx = nx
        self.ny = ny
        
        # 双缓冲
        self.E_current = np.zeros((nx, ny))
        self.E_next = np.zeros((nx, ny))
        self.H_current = np.zeros((nx, ny))
        self.H_next = np.zeros((nx, ny))
    
    def step_overlapped(self):
        # 1. 计算内部点（不依赖边界）
        self.compute_interior(self.E_current, self.H_current, 
                             self.E_next, self.H_next)
        
        # 2. 启动非阻塞边界交换
        requests = self.start_boundary_exchange(self.E_current)
        
        # 3. 计算边界点（与通信重叠）
        self.compute_boundary(self.E_current, self.H_current,
                             self.E_next, self.H_next)
        
        # 4. 等待通信完成
        MPI.Request.Waitall(requests)
        
        # 5. 交换缓冲区
        self.E_current, self.E_next = self.E_next, self.E_current
        self.H_current, self.H_next = self.H_next, self.H_current

工程应用案例

案例1：飞机RCS计算的并行加速

　　问题描述：
计算F-16战斗机在X波段（10 GHz）的雷达散射截面，目标尺寸约15m，需要约10亿个网格单元。

　　并行策略：

class AircraftRCS_Simulation:
    def __init__(self, aircraft_model, frequency, n_processes):
        self.frequency = frequency
        self.wavelength = 3e8 / frequency
        self.n_processes = n_processes
        
        # 网格划分：λ/10分辨率
        dx = self.wavelength / 10
        self.nx = int(aircraft_model.length / dx)
        self.ny = int(aircraft_model.wingspan / dx)
        self.nz = int(aircraft_model.height / dx)
        
        print(f"Grid size: {self.nx} x {self.ny} x {self.nz}")
        print(f"Total cells: {self.nx * self.ny * self.nz / 1e9:.2f} billion")
    
    def setup_parallel(self):
        """
        设置并行计算环境
        """
        # 使用3D域分解
        px = py = pz = int(round(self.n_processes ** (1/3)))
        
        # 创建笛卡尔拓扑
        self.fdtd = FDTD3D_MPI(self.nx, self.ny, self.nz, 
                               n_steps=5000, px=px, py=py, pz=pz)
    
    def compute_rcs(self, theta_range, phi_range):
        """
        计算双站RCS
        """
        rcs_results = {}
        
        for theta in theta_range:
            for phi in phi_range:
                # 设置入射波方向
                self.set_incident_wave(theta, phi)
                
                # 运行FDTD
                self.fdtd.run()
                
                # 近场到远场变换
                far_field = self.near_to_far_field_transform()
                
                # 计算RCS
                rcs = 4 * np.pi * np.abs(far_field)**2
                rcs_results[(theta, phi)] = rcs
        
        return rcs_results
    
    def benchmark(self):
        """
        性能基准测试
        """
        import time
        
        times = []
        for n_procs in [1, 4, 16, 64, 256]:
            start = time.time()
            # 运行简化测试
            elapsed = time.time() - start
            times.append((n_procs, elapsed))
            
            efficiency = times[0][1] / (n_procs * elapsed)
            print(f"Processes: {n_procs}, Time: {elapsed:.2f}s, Efficiency: {efficiency:.2%}")

案例2：5G城市传播的大规模射线追踪

　　问题描述：
在1km²城市环境中，使用射线追踪预测5G毫米波（28 GHz）覆盖，包含数千栋建筑物。

　　并行实现：

class ParallelRayTracing:
    def __init__(self, city_model, frequency, n_rays=1e6):
        self.city = city_model
        self.frequency = frequency
        self.n_rays = int(n_rays)
        self.wavelength = 3e8 / frequency
    
    def parallel_trace(self, tx_position, n_processes):
        """
        并行射线追踪
        """
        from mpi4py import MPI
        
        comm = MPI.COMM_WORLD
        rank = comm.Get_rank()
        size = comm.Get_size()
        
        # 每个进程处理一部分射线
        rays_per_process = self.n_rays // size
        start_ray = rank * rays_per_process
        end_ray = start_ray + rays_per_process
        
        local_results = []
        
        for ray_id in range(start_ray, end_ray):
            # 生成射线
            direction = self.generate_ray_direction(ray_id)
            
            # 追踪射线
            path = self.trace_ray(tx_position, direction)
            
            # 计算接收功率
            received_power = self.compute_path_loss(path)
            
            local_results.append({
                'ray_id': ray_id,
                'path': path,
                'power': received_power
            })
        
        # 收集所有结果
        all_results = comm.gather(local_results, root=0)
        
        if rank == 0:
            # 合并结果
            coverage_map = self.build_coverage_map(all_results)
            return coverage_map
    
    def gpu_accelerated_intersection(self, rays, buildings):
        """
        GPU加速的射线-建筑物相交检测
        """
        import cupy as cp
        
        # 将数据转移到GPU
        rays_gpu = cp.array(rays)
        buildings_gpu = cp.array(buildings)
        
        # CUDA核函数进行相交检测
        intersections = self.cuda_intersection_kernel(rays_gpu, buildings_gpu)
        
        return cp.asnumpy(intersections)
    
    def hybrid_computation(self):
        """
        CPU+GPU混合计算
        """
        # GPU处理主要射线追踪
        gpu_rays = int(self.n_rays * 0.8)
        gpu_results = self.gpu_ray_tracing(gpu_rays)
        
        # CPU处理复杂反射/绕射
        cpu_rays = self.n_rays - gpu_rays
        with multiprocessing.Pool() as pool:
            cpu_results = pool.map(self.cpu_trace_worker, range(cpu_rays))
        
        return self.merge_results(gpu_results, cpu_results)

案例3：超材料设计的参数扫描优化

　　问题描述：
优化超材料单元结构，需要在参数空间进行大规模扫描。

　　并行优化：

class ParallelOptimization:
    def __init__(self, design_space, objective_func):
        self.design_space = design_space
        self.objective = objective_func
    
    def grid_search_parallel(self, n_processes):
        """
        并行网格搜索
        """
        from mpi4py import MPI
        import itertools
        
        comm = MPI.COMM_WORLD
        rank = comm.Get_rank()
        size = comm.Get_size()
        
        # 生成参数组合
        param_grid = list(itertools.product(*self.design_space))
        
        # 分配任务
        local_params = param_grid[rank::size]
        
        local_results = []
        for params in local_params:
            # 评估目标函数
            result = self.objective(params)
            local_results.append((params, result))
        
        # 收集结果
        all_results = comm.gather(local_results, root=0)
        
        if rank == 0:
            # 找到最优解
            best_result = min([r for sublist in all_results for r in sublist], 
                            key=lambda x: x[1])
            return best_result
    
    def genetic_algorithm_parallel(self, population_size, generations):
        """
        并行遗传算法
        """
        from mpi4py import MPI
        import numpy as np
        
        comm = MPI.COMM_WORLD
        rank = comm.Get_rank()
        size = comm.Get_size()
        
        # 每个进程维护一部分种群
        local_pop_size = population_size // size
        
        for gen in range(generations):
            # 本地适应度评估
            local_fitness = self.evaluate_population(local_population)
            
            # 全局最优共享
            global_best = comm.allreduce(local_best, op=MPI.MIN)
            
            # 选择、交叉、变异
            local_population = self.evolve_population(local_population, 
                                                      local_fitness,
                                                      global_best)
            
            # 迁移（定期交换个体）
            if gen % 10 == 0:
                self.migrate_individuals(local_population)
        
        return global_best

总结与展望

关键知识点总结

1. 并行计算基础：
   • 理解并行计算模型（共享内存vs分布式内存）
   • 掌握Amdahl定律与可扩展性分析
   • 熟悉现代处理器架构（多核CPU、GPU）
2. OpenMP并行化：
   • 循环级并行与任务并行
   • 线程同步与数据竞争避免
   • SIMD向量化优化
3. MPI分布式计算：
   • 域分解与负载均衡
   • 点对点与集合通信
   • 通信优化与计算重叠
4. GPU加速：
   • CUDA编程模型与内存管理
   • 核函数优化与线程配置
   • 多GPU并行与NCCL通信
5. 性能优化：
   • 内存层次优化（缓存、分块）
   • 计算与通信重叠
   • 异构计算策略

发展趋势

1. E级计算（Exascale）：
   • 10^18 FLOPS计算能力
   • 百亿亿次网格处理能力
   • 新型互连网络与存储
2. AI加速器集成：
   • TPU、NPU等专用加速器
   • 神经网络代理模型
   • 物理信息神经网络（PINN）
3. 量子计算探索：
   • 量子电磁学算法
   • 量子-经典混合计算
   • 量子优势应用识别

参考文献

1. Gropp, W., Lusk, E., & Skjellum, A. (1999). Using MPI. MIT Press.
2. Chapman, B., Jost, G., & Van Der Pas, R. (2007). Using OpenMP. MIT Press.
3. Sanders, J., & Kandrot, E. (2010). CUDA by Example. Addison-Wesley.
4. Taflove, A., & Hagness, S. C. (2005). Computational Electrodynamics. Artech House.
5. Yu, W., et al. (2006). Parallel Finite-Difference Time-Domain Method. Artech House.

附录：性能优化检查清单

优化级别	检查项	预期收益
算法	选择合适算法复杂度	10-1000x
实现	循环优化、向量化	2-10x
内存	缓存优化、数据对齐	2-5x
并行	OpenMP/MPI并行化	接近核心数
加速	GPU加速	10-100x
通信	重叠、聚合	1.5-3x