简介面向光子逆设计与拓扑优化研究者的 Python 基准测试套件基于 ceviche 有限差分频域模拟器与 HIPS autograd 自动微分构建内置波导分束器、模式转换器、波导弯曲和波分复用器等典型集成光子组件可对不同拓扑优化算法进行横向对比与性能评估。压缩包共 43 个文件以 35 个 Python 脚本为主涵盖散射参数计算、模式与基础原语定义、模型封装及测试用例另附 4 张组件示意图与 README 文档整体约 772KB。该套件源自已发表论文的配套代码可复现性强目前已有 369 人学习下载。借助现有测试脚本可快速搭建环境模块化目录结构便于替换优化策略或调整组件参数开展自定义逆向设计实验下载后从 README.md 即可掌握整体组织与调用方式。1. 光子逆设计基准测试为什么你的拓扑优化算法需要一个公共赛道有一次我用自己写的拓扑优化程序做波导模式转换器在某个测试算例上收敛得又快又稳换了一组边界条件就整个发散。排查到最后发现问题根本不在算法而在仿真设置有的挑战给的是点源有的是模式激发网格分辨率也不统一。这时候我强烈感觉到算法之间缺少一个公共“比赛场地”。ceviche-challenges 就是这套 Python 代码的定位把分束器、模式转换器、波导弯曲、波分复用器这些经典集成光子器件封装成标准逆设计挑战所有挑战共享同一套有限差分频域FDFD模拟接口并用 autograd 把梯度暴露给优化器。无论你是做拓扑优化算法的人还是刚进入光子逆向设计的开发者都可以把它当成公共赛道用同一把尺子量不同算法的真实水平。2. 技术栈与工程结构FDFD模拟器、自动微分和挑战套件的配合2.1 为什么选FDFD当仿真底座光子逆设计本质上是把仿真变成可优化的计算图。频域仿真方法有 FDFD、FDTD、本征模展开几种FDFD 的特别之处在于它把 Maxwell 方程离散成稀疏线性系统求解一次频域响应不需要时间步进。对设计参数求梯度时也能直接对线性系统做灵敏度分析配合自动微分可以省掉手写伴随求解器的工作量。ceviche 是一个很轻量的 FDFD 实现yee 网格离散、PML 吸收边界、模式激发都封装好了特别适合作为逆向设计的底层。用 ceviche 当底座还有一个实际好处它的核心操作都是 numpy 原生的能直接喂给 HIPS autograd。常见做法是把整个仿真流程写成一个大的 numpy 函数输入是设计密度输出是散射参数最后通过一个标量损失驱动梯度下降。因为 FDFD 不需要逐步迭代时间计算图相对紧凑反向传播的中间变量也不像 FDTD 那样占内存。实际跑起来一个普通设计域在单卡 GPU 上做一二百轮优化是完全可接受的。我比较看重的是它把“可微仿真”做成了标准 API而不是散落的脚本。工程上这意味着新算法可以直接调用统一的 loss 和 grad 接口不需要每个研究者都重写底层求解器。再加上 autograd 对 numpy 子集的支持梯度计算在大多数情况下是“白盒”的出错了也可以从执行轨迹里排查而不是把仿真器当成黑匣子。2.2 目录结构挑战套件如何分层解压下载的压缩包后根目录是 ceviche-challenges-main。核心包是 ceviche_challenges里面文件按作用可以分成四层。第一层是仿真底座包括 scattering.py、modes.py、ops.py、units.py第二层是挑战定义包括 model_base.py、params.py、primitives.py、defs.py第三层是四个具体器件模块beam_splitter、mode_converter、waveguide_bend、wdm第四层是测试文件scattering_test.py、modes_test.py、ops_test.py、primitives_test.py、defs_test.py 全部位于 ceviche_challenges 包内用 pytest 可以直接运行。目录结构如下ceviche-challenges-main/ ├── ceviche_challenges/ │ ├── model_base.py │ ├── params.py │ ├── primitives.py │ ├── ops.py │ ├── units.py │ ├── scattering.py │ ├── modes.py │ ├── defs.py │ ├── beam_splitter/ │ ├── mode_converter/ │ ├── waveguide_bend/ │ ├── wdm/ │ ├── scattering_test.py │ ├── modes_test.py │ ├── ops_test.py │ ├── primitives_test.py │ └── defs_test.py ├── img/ │ ├── beam_splitter.png │ ├── mode_converter.png │ ├── waveguide_bend.png │ └── wavelength_divison_multiplexer.png ├── setup.py └── README.md这种分层的思路很清晰。仿真底座不动挑战定义层把器件目标描述成参数空间里的损失函数测试文件同时承担“环境自检”和“数值正确性验证”两种角色。对二次开发来说新器件只需要在 defs.py 或新增子模块里实现一个挑战类不需要动仿真底座。setup.py 负责包安装README.md 里会写明具体用法img 目录里的 png 是每个挑战问题的目标版图示例。2.3 参数体系、单位制与制造约束逆设计通常把设计区域离散成 Nx×Ny 的像素网格每个像素对应一个密度值取值在 0 到 1 之间。0 代表背景材料通常是二氧化硅1 代表核心材料比如硅。params.py 负责维护这个密度数组并对外暴露形状、边界条件、初始值等属性。modes.py 则负责把端口定义成具体的波导模式这样仿真时就知道从哪一侧输入、从哪一侧接收。units.py 里定义了单位换算。FDFD 是在物理尺度上求解因此波长、波导宽度、设计区域尺寸都必须统一。比如模式转换器的设计域常见尺寸是 2.4μm × 0.9μm工作波长 1550nm。这套代码跑出的结果可以直接映射到光刻版图这是我觉得比较省心的一点。很多教学代码用的都是无量纲单位真到流片环节就翻车还得重新缩放一堆参数。制造约束体现在 ops.py 里的投影操作。裸密度优化会出现大量 0.5 左右的中间值代工厂没法刻写。常见做法是先对密度做高斯模糊再用双曲正切投影把密度推向 0 或 1。这个投影不是一次性完成的通常在优化过程中让 beta 参数从 1 慢慢升到几十早期给算法足够探索空间后期强制执行二值化。下面这段代码是投影函数最典型的写法import numpy as np def projection(density, eta0.5, beta1.0): # 将密度投影到接近0或1的二值分布 # eta: 阈值通常取0.5 # beta: 投影强度越大越接近二值 numerator np.tanh(beta * eta) - np.tanh(beta * (density - eta)) denominator np.tanh(beta * eta) - np.tanh(beta * (1 - eta)) return numerator / denominator逻辑说明density 是原始设计变量取值范围 [0,1]eta 是阈值默认 0.5beta 控制投影强度。当 beta 很小时这个函数近似恒等映射保留连续密度当 beta 增大时小于 eta 的值被压向 0大于 eta 的值被压向 1。需要注意 beta 太大会导致梯度消失所以我一般会做 beta 退火每跑几十轮就把 beta 乘以 1.1而不是一开始就设成 50。参数解释density 的形状就是设计域的像素网格形状eta 不要随意改成接近 0 或 1 的值否则投影会严重偏置优化结果容易陷入局部极小。3. 核心实现拆解一个挑战问题从定义到梯度下降的完整链路3.1 model_base.py统一挑战接口整个套件能称为“基准测试”靠的是 model_base.py 里定义的统一基类。每个具体器件子类都继承这个基类并实现 forward、loss、target 这几个方法。设计者只需要把注意力放在自己的算法上挑战问题的物理定义全部封装在类内部。这里我按最常见的写法补一段基类骨架import numpy as np from autograd import grad class ChallengeBase: def __init__(self, params, wavelengths, portsNone): self.params params # 设计变量容器包含网格形状 self.wavelengths wavelengths # 工作波长数组 self.ports ports or [] # 端口定义来自 modes.py def forward(self, design): # 输入一维设计密度返回S参数数组 raise NotImplementedError def target(self): # 返回目标S参数子类必须实现 raise NotImplementedError def loss(self, design): # 均方误差损失 s_actual self.forward(design) return np.mean((s_actual - self.target()) ** 2) def gradient(self, design): # 用 autograd 对 loss 求梯度 return grad(self.loss)(design)逻辑说明forward 是仿真链路入口内部会调用 ops 把密度转成介电常数再调用 scattering 仿真得到 S 参数。loss 直接算仿真结果和目标响应的均方误差gradient 用 autograd 对 loss 自动微分。这个设计的巧妙之处在于子类只需要实现 forward 和 target梯度是自动拿到的不会因为每个器件目标不同而重复写伴随求解器。参数说明design 是长度为 Nx*Ny 的展平数组因为 autograd 对高维数组的索引操作有时会拖慢速度展平后更稳wavelengths 可以是标量也可以是数组如果做 WDM 这类多波长器件通常要传一组波长loss 里会对多个波长求平均。3.2 four器件分束器、模式转换器、波导弯曲与WDM四个器件分别放在 beam_splitter、mode_converter、waveguide_bend、wdm 子目录里。它们的目标各不相同但都通过继承 ChallengeBase 来统一接口。以 mode_converter 为例它的目标是把输入波导的基模 TE0 模式转换到输出波导的 TE1 模式传统波导没法直接在短距离内完成这种转换拓扑优化可以在设计域里挖出任意形状的结构。模式转换器的挑战定义核心是 target() 方法。常见实现是这样的class ModeConverter(ChallengeBase): def __init__(self): # 设计域尺寸 2.4um x 0.9um像素网格 120x45 params Params(shape(120, 45), physical_size(2.4, 0.9)) super().__init__(paramsparams, wavelengths[1.55]) def target(self): # 单波长目标传输端 S21 到 0.95反射端 S11 到 0.05 return np.array([0.95, 0.05]) def forward(self, design): epsr ops.density_to_epsr(design, background1.0, core12.25) return scattering.simulate_s_params(epsr, self.wavelengths[0], self.mode_ports)参数说明params 里的 shape 是像素网格的行列数physical_size 是微米级的物理尺寸wavelengths 单位也是微米1.55 就是 1550nm。target 中的 S21 是传输系数S11 是反射系数目标值不是严格 1.0是为了给工艺误差留余量。beam_splitter 的目标则是让输出两个端口的功率各占一半S21 和 S31 都接近 0.707waveguide_bend 要在紧凑转角内维持低反射wdm 要按波长把光分开到不同端口所以它的 target 是多波长数组例如波长 λ1 时端口 A 高透射、端口 B 低透射波长 λ2 时相反。3.3 ops.py与primitives.py可微几何生成拓扑优化最容易被忽略的是“密度到物理结构”的映射。设计变量只是个 0 到 1 的数组但仿真需要介电常数分布而且代工厂需要结构不悬空、不断裂。ops.py 里除了上一章提到的 projection还有 blur 和高斯滤波。blur 的作用是让相邻像素的密度平滑过渡避免优化出的结构出现孤立点。primitives.py 是生成初始结构的地方。比如波导弯曲挑战输入和输出波导的位置是固定的设计域中间是空白算法需要自己生成一个连接两个端口的弯曲形状。primitives 提供直线、矩形、贝塞尔曲线的可微版本这样初始设计也可以参与优化而不是一开始就全零。下面是一个可微模糊滤波的实现思路from scipy.ndimage import gaussian_filter from autograd.extend import primitive, defvjp primitive def blur(density, sigma): # 对密度做高斯模糊sigma控制平滑半径 return gaussian_filter(density, sigmasigma) def blur_vjp(ans, density, sigma): # 高斯卷积是对称算子VJP就是自己对自己 def vjp(g): return gaussian_filter(g, sigmasigma), None return vjp defvjp(blur, blur_vjp)逻辑说明这里把 scipy 的 gaussian_filter 包装成 autograd 的 primitive并定义了它的向量-雅可比积。因为高斯卷积是自伴随的反向传播时只需要对上游梯度再卷一次所以 vjp 函数里直接调用 gaussian_filter。参数说明sigma 是模糊半径单位是像素太小起不到连接作用太大又会让结构边界模糊我一般取 1.5 到 2.5。注意如果直接用 scipy 的 gaussian_filter 而不注册 vjpautograd 会在反向传播时报错这也是很多新手第一次跑这个套件时最常见的坑。3.4 测试文件数值梯度验证与仿真正确性测试文件是这套代码里被低估的部分。ops_test.py 会验证 blur 和 projection 的数值梯度与 autograd 梯度是否一致方法是用有限差分法算梯度然后和自动微分结果对比。scattering_test.py 会自动创建一个直波导检查 S21 是否接近 1这一条能一次性验证 FDFD 模拟器、模式激发、端口数据处理整条链路是否正常。defs_test.py 则会检查每个挑战问题的 target 参数是否在物理合理范围内比如分束器的两个输出端口功率和是否接近 1。这些测试不是摆设做逆向设计时最怕底层仿真与梯度不自洽。数值梯度校验如果通不过后面再花哨的优化算法都跑不出可信结果。我自己的习惯是拿到这份代码后第一件事就是跑 pytest全绿后才会开始调算法。4. 运行实战与避坑从安装到误差控制4.1 环境准备用虚拟环境隔离依赖这套 Python 代码依赖 numpy、scipy、autograd、ceviche 几个包直接 pip install -e . 会把项目本身安装到当前环境。但 ceviche 在部分 Python 版本下存在编译兼容问题所以我建议创建独立的虚拟环境避免污染其他项目。安装步骤如下python -m venv ceviche-env source ceviche-env/bin/activate pip install --upgrade pip pip install -e .如果执行 pip install -e . 后报找不到 ceviche需要手动补装pip install ceviche。装完后再运行 python -c import ceviche 确认导入成功。要注意 autograd 与 numpy 的版本组合偶尔会出现 API 兼容问题常见现象是 import 时 AttributeError解决方法是把 numpy 降到较旧的稳定版本比如 1.20 左右再重新安装 autograd。这个组合对 Python 3.7 到 3.9 都算比较稳。4.2 跑通基准测试从pytest到自定义优化脚本环境准备好之后先跑测试套件确认所有组件正常pip install pytest pytest ceviche_challenges -v这里会执行 scattering_test、ops_test、primitives_test、defs_test 等所有测试。如果全部通过说明 FDFD 仿真和自动微分链路是自洽的。接下来我通常写一个最小优化脚本验证一个挑战问题的 loss 能不能被梯度下降拉低import numpy as np from autograd import grad from ceviche_challenges import mode_converter model mode_converter.Challenge() design np.random.uniform(0.4, 0.6, sizemodel.params.shape) loss_fn model.loss grad_fn grad(loss_fn) for i in range(100): g grad_fn(design) design - 0.05 * g design np.clip(design, 0, 1) if i % 10 0: print(fiter {i}, loss {loss_fn(design):.6f})参数说明design 初始值取 0.4 到 0.6 之间的随机数而不是 0 到 1 均匀随机这样仿真一开始就有一个介电常数不为零的物理结构避免 FDFD 矩阵奇异。learning rate 0.05 只是经验值实际要根据损失尺度和梯度范数调整。每一步对设计变量做 clip 是为了满足制造约束。注意这里不能直接用 model.params 里的数据因为 params 可能被内部类用于缓存需要单独创建一个 ndarray 传给 loss。4.3 常见问题排查依赖冲突、梯度NaN与端口定义问题一import 阶段报 AttributeError: module numpy has no attribute float_。现象是 pip install 完成后导入 autograd 或 ceviche 时报错。原因是最新 numpy 移除了 numpy.float_ 别名而旧版 ceviche 或 autograd 内部还在引用它。解决方法是把 numpy 降级到 1.23 或更低或者升级到能兼容的 ceviche 版本。这类问题没有黑匣子可看直接看 traceback 最后一行就能定位到是哪个包在调用。问题二ops_test.py 数值梯度失败误差在 1e-2 数量级。现象是测试打印有限差分梯度与 autograd 梯度不一致。原因通常是自定义 primitive 的 vjp 写错了或者模糊半径 sigma 随机初始化导致高斯近似不准确。解决方法是先固定 sigma 为一个较小整数比如 1.0再做梯度对比如果误差仍大检查 blur_vjp 里是否漏掉了对 sigma 的梯度虽然 sigma 一般不需要梯度但 primitive 定义里切不可省略参数占位。问题三loss 在优化初期就变成 NaN。现象是第 1 轮迭代 loss 正常第 5 轮突然全 NaN。原因一般是设计变量某处变成极端值导致介电常数分布出现零或负值FDFD 主对角矩阵奇异。解决方法是给投影函数传入一个安全的 epsilon并在每次迭代后检查 density 的 min 和 max如果 min 0 或 max 1用 clip 强制拉回并降低学习率。这个问题我遇到过不止一次最后发现是 blur 在边界处的 padding 方式导致的需要用 modenearest 而不是默认的 modereflect。问题四不同机器上跑出的拓扑结构差异很大。现象是同样的参数两台机器跑出的 loss 曲线相近但最终版图完全不同。原因有两个一是 numpy 随机种子没有固定二是 FDFD 的稀疏求解器用到了 OpenBLAS 多线程线程数不同导致浮点累加顺序变化。解决方法是脚本开头设置 np.random.seed(0)同时设置环境变量 OPENBLAS_NUM_THREADS1 和 OMP_NUM_THREADS1让求解可复现。这类数值稳定性问题在拓扑优化里很常见不是算法错了是环境差异。问题五S21 大于 1物理上不合理。现象是某个端口透射系数出现 1.2 之类的值。原因多半是端口模式功率归一化没有做或者模式激发时没有减去反射功率。解决方法是查看 modes.py 里模式投影的写法确认使用的是归一化场分布而不是原始场强。这一条不是套件本身的 bug而是我改造端口时踩过的坑如果你自己动了 modes 相关的代码一定跑一遍 scattering_test 验证直波导透射约等于 1。5. 二次开发把一个新的光电器件加入挑战套件如果你想加入一个现有四个器件之外的挑战比如 Y 型分束器关键不是写仿真器而是写一个继承 ChallengeBase 的子类。先把 Y 型分束器的物理目标量化输入端口接 TE0 模式输出端是两个对称分支目标是在 1550nm 处两个输出端口的功率各占一半反射系数小于 0.1。在代码里这就是一个 target 数组为 [0.71, 0.71, 0.1] 的类。先创建子包 y_branch在init.py 里写挑战类import numpy as np from ..model_base import ChallengeBase from .. import ops, scattering, modes class YBranchChallenge(ChallengeBase): def __init__(self): params Params(shape(120, 120), physical_size(2.4, 2.4)) super().__init__(params, wavelengths[1.55]) def target(self): return np.array([0.71, 0.71, 0.1]) def forward(self, design): epsr ops.density_to_epsr(design, background1.0, core12.25) s scattering.simulate_s_params(epsr, self.wavelengths[0], ports) # 假设返回 [S21, S31, S11] return np.array([s[1], s[2], s[0]])然后写一个测试文件 test_y_branch.py内容先验证梯度一致性再验证随机初始设计的 loss 不为 NaN。这样做能保证新挑战的仿真链路和梯度链路都是通的不会污染原有套件。从那以后我每次加新挑战都会先写一个最小优化脚本把初始设计设成常量跑 20 步梯度下降确认 loss 是在下降而不是抖动然后才去调材料和边界参数。否则你很可能在梯度本身就缺失的情况下花一整天调优化器参数最后还以为是算法出了问题。这套路和调卷积网络一样先跑通一条小曲线再谈精度。希望这个思路能帮你在自己的器件设计里少走几步弯路。本文还有配套的精品资源点击获取
