meent开源库实战：RCWA/TMM原理、实现与超表面优化避坑指南

1. 项目概述与核心价值如果你正在设计光子晶体、超表面或者任何带有周期性微纳结构的光学器件那么“仿真”这一步几乎是绕不开的。无论是想优化一个光栅耦合器的耦合效率还是设计一个能将特定波长光高效偏转的衍射元件你都需要一个可靠的工具来预测光与结构相互作用的结果。在众多电磁仿真方法中严格耦合波分析RCWA因其在处理周期性结构时的天然优势和计算效率成为了这个领域的“标准答案”之一。然而从教科书上的矩阵方程到真正能跑出可靠结果的代码中间隔着巨大的鸿沟——数值稳定性、计算效率、以及如何将复杂的数学公式转化为清晰可维护的程序逻辑这些都是实践中会遇到的真实挑战。最近我在一个超表面优化项目中深入使用了meent这个开源Python库。它完整实现了RCWA结合传输矩阵法TMM的求解流程并且原生支持自动微分AD这对于基于梯度的器件逆向设计来说简直是“神器”。更重要的是它的代码结构相当清晰把卷积矩阵生成、特征值求解、场连接等核心步骤都模块化了非常适合用来学习和理解RCWA的底层实现。本文将结合我的使用经验带你从RCWA/TMM的核心原理出发一步步拆解meent的实现并分享在实战中关于精度、速度和稳定性的那些“坑”与技巧。无论你是刚接触计算电磁学的学生还是需要快速上手一个可靠仿真工具的研究者相信都能从中找到有价值的内容。2. RCWA与TMM从物理方程到矩阵求解2.1 核心思想在傅里叶空间中求解麦克斯韦方程RCWA的基本思想非常直观对于一个在x和y方向具有周期性的结构光栅其中的电磁场也可以展开成具有相同空间频率的傅里叶级数。这就是弗洛奎特-布洛赫定理的体现。我们将介电常数分布ε(x, y)也进行傅里叶展开。这样时谐麦克斯韦方程组中的旋度方程在傅里叶空间中就从微分方程变成了代数方程。具体来说我们从麦克斯韦旋度方程出发∇ × E -jωμH和∇ × H jωεE。 在傅里叶空间中空间微分算符∇变成了波矢k的乘法运算。对于周期性结构k可以写为入射波矢k_inc加上倒格矢的整数倍(n, m)即k_x, k_y构成对角矩阵Kx和Ky。介电常数及其倒数在傅里叶空间中的表示则是一个托普利茨Toeplitz矩阵也称为卷积矩阵记作[ε]和[1/ε]。将电场E和磁场H的傅里叶级数系数向量代入方程经过一系列推导核心是消去z方向的导数得到关于横向场分量的方程最终可以化归为一个标准的特征值问题∂^2/∂z^2 [E_x, E_y]^T -k0^2 * A * [E_x, E_y]^T其中矩阵A由Kx,Ky,[ε],[1/ε]等构成。求解这个特征值问题我们就能得到每一层光栅中可能存在的传播模式特征模及其在z方向的传播常数特征值。注意符号约定的重要性在推导中一个非常容易出错的地方是符号约定。例如原文公式中提到的负号调整是为了确保旋度运算E × H的方向与波传播方向k_g,z一致。不同的文献可能采用e^(jωt)或e^(-iωt)的时谐约定这会影响到所有虚数单位j的符号进而影响特征值问题的矩阵构造。meent内部采用e^(-iωt)约定这是物理学中的常见约定如果你从采用e^(jωt)约定的工程文献推导公式直接套用可能会导致结果错误。在阅读代码和公式时务必先确认其时空谐约定。2.2 层间连接传输矩阵法TMM的精髓求解出单层的特征模后我们得到了该层内上行波和下行波的系数。但一个器件通常是多层结构比如基底、光栅层、覆盖层。TMM的作用就是优雅地处理这些层之间的边界条件。在每一层的界面处电磁场的切向分量E_x, E_y, H_x, H_y必须连续。我们将这些条件用矩阵形式写出。对于第ℓ层其界面上的场可以用该层的特征向量矩阵W_ℓ、特征值矩阵q_ℓ以及层厚度d_ℓ构成的相位矩阵X_ℓ来表示。以原文中的公式为例在输入边界z0和输出边界zd处我们可以建立联系入射场I、反射场系数R、透射场系数T以及各层内部模式系数c^±的大型线性方程组。通过消去内部系数c^±最终可以得到一个简洁的传输矩阵关系[I; R] M_total * [T; 0]其中M_total是所有层传输矩阵的连乘。从这个关系式中我们就能直接解出我们关心的反射系数R和透射系数T。2.3 衍射效率的计算从复振幅到真实功率得到复振幅系数R和T后我们最终需要的是可测量的物理量——衍射效率DE即衍射到某一级次的光功率与入射光功率的比值。这里需要用到坡印廷矢量进行时间平均计算。对于第(n, m)级衍射其反射效率DE_r和透射效率DE_t的计算公式如原文所示。公式中包含了Re(k_z / (k0 n cosθ))这样的因子它本质上是z方向波矢的归一化确保了对于倏逝波k_z为虚数其衍射效率为零因为倏逝波不携带能量传播。这是RCWA计算中一个非常重要的自洽性检查点所有传播级次k_z为实数的衍射效率之和应等于1无损耗时或小于1有材料吸收时。3. meent库实战从安装到第一个仿真3.1 环境搭建与基础概念meent可以通过pip直接安装pip install meent。它支持NumPy、JAX和PyTorch三种后端。选择哪个后端取决于你的需求NumPy (backend0)最稳定兼容性最好但不支持自动微分AD。JAX (backend1)支持AD并且通过JIT即时编译可以获得不错的性能但在处理超大矩阵高傅里叶级数时首次编译开销和内存管理可能成为瓶颈。PyTorch (backend2)支持AD生态完善GPU支持良好是进行梯度优化研究的推荐选择。初始化一个仿真实例非常简单import meent import numpy as np # 定义一个一维光栅结构2层每层10个像素 ucell np.array([ [[1, 1, 1, 3.48, 3.48, 3.48, 3.48, 1, 1, 1]], # 第一层空气(1)和硅(3.48) [[1, 3.48, 3.48, 1, 1, 1, 1, 3.48, 3.48, 1]], # 第二层 ]) # 形状为 (层数, Y方向像素数, X方向像素数)。对于1D光栅Y方向为1。 # 创建仿真器 mee meent.call_mee( backend0, # 使用NumPy后端 grating_type0, # 1D光栅非锥形入射 pol0, # TE偏振 (pol0 - ψ90°) n_I1.0, # 上覆层通常是空气折射率 n_II1.45, # 衬底例如二氧化硅折射率 theta0., # 入射角弧度0表示正入射 phi0., # 方位角弧度 wavelength900, # 真空波长单位nm fourier_order10, # 傅里叶截断阶数实际使用21个谐波-10到10 period[1000], # 光栅周期单位nm ucellucell, # 定义的结构 thickness[200, 300], # 每层的厚度单位nm顺序从上到下 fft_type0, # 使用离散傅里叶级数(DFS)进行栅格建模 )这段代码定义了一个简单的双层一维光栅。ucell是一个三维数组直接定义了每个像素点的介电常数或折射率的平方。fft_type0表示使用离散傅里叶变换DFS来处理这种像素化的结构。3.2 两种结构建模方式栅格与矢量meent提供了两种定义结构的方式对应不同的应用场景和精度需求。栅格建模 (Raster Modeling,fft_type0 或 1)这就是上面例子中使用的方式。你将设计区域离散化为一个个小像素体素并为每个像素赋值介电常数。这种方式非常灵活可以描述任意形状是拓扑优化Topological Optimization的天然选择因为每个像素都可以作为一个独立的设计变量。fft_type0: 使用离散傅里叶级数DFS。计算速度快但当像素尺寸与特征尺寸相比不够精细时会产生“阶梯误差”和吉布斯现象导致精度下降。fft_type1: 使用连续傅里叶级数CFS处理栅格。它通过解析积分计算傅里叶系数精度高于DFS尤其对于高对比度结构但计算量稍大且当前版本不支持反向传播AD。矢量建模 (Vector Modeling,fft_type2)这种方式直接描述几何形状的边界。例如你可以定义一个旋转矩形、椭圆等。meent内部会将这些形状分解为多个与坐标轴对齐的子矩形然后利用CFS和sinc函数精确计算其傅里叶系数。import torch # 使用PyTorch后端以支持后续梯度计算 mee meent.call_mee(backend2, grating_type2, fft_type2, period[1000, 1000], ...) # 定义一个旋转矩形作为结构 length_x torch.tensor(400.0, requires_gradTrue) # 将尺寸设为可求导变量 length_y torch.tensor(600.0, requires_gradTrue) angle torch.tensor(30 * np.pi / 180, requires_gradTrue) center [500, 500] # 单元中心坐标 n_split [20, 20] # 将形状近似为20x20个小矩形精度更高 n_index 3.48 # 结构材料折射率硅 # 创建旋转矩形对象 rect_obj mee.rectangle_rotate(*center, length_x, length_y, *n_split, n_index, angle) # 定义层信息基底折射率 对象列表 layer_info_list [[1.0, [rect_obj]]] # 基底为空气(1.0)上面有一个硅矩形 mee.draw(layer_info_list) # 可视化结构矢量建模的精度通常高于栅格建模特别是对于光滑边界。它适用于形状导数Shape Derivative优化即优化结构的尺寸、位置、旋转角度等参数而保持拓扑不变。n_split参数控制着近似精度值越大分解的子矩形越多形状越精确但计算量也越大。实操心得如何选择建模方式做拓扑优化像素级变化无脑选栅格建模 (fft_type0)。你可以方便地将ucell的每个元素作为优化变量。做形状/尺寸优化如优化圆柱半径、光栅占空比优先选矢量建模 (fft_type2)。它得到的梯度形状导数更准确且设计变量更少优化问题更简洁。追求最高精度进行验证计算如果结构是简单几何形状用矢量建模高n_split。如果是复杂形状可以尝试栅格建模fft_type1(CFS)并与商业软件结果交叉验证。处理结构重叠在矢量建模中layer_info_list中对象的顺序决定了堆叠层次。列表后面的对象会覆盖前面的对象。这为你设计复杂多层结构提供了灵活性。3.3 运行仿真与结果提取设置好结构后运行仿真只需要一行代码de_ri, de_ti mee.conv_solve()de_ri和de_ti分别是反射和透射的衍射效率数组。对于1D光栅它们是一维数组索引对应衍射级次。例如de_ti[10]对应透射的10级衍射效率。0级通常对应镜面反射或直接透射。如果你想进一步分析器件内部的电磁场分布可以使用conv_solve_field方法或先conv_solve再calculate_field# 方法1一次性计算效率和场 de_ri, de_ti, field_cell mee.conv_solve_field(res_x50, res_y50, res_z30) # 方法2分两步 de_ri, de_ti mee.conv_solve() field_cell mee.calculate_field(res_x50, res_y50, res_z30)field_cell是一个四维数组其形状和含义取决于光栅类型1D光栅 (TE或TM)形状为(总Z层数, res_y, res_x, 3)。最后一维的3个分量分别是TE偏振:[Ey, Hx, Hz]TM偏振:[Hy, Ex, Ez]1D锥形入射或2D光栅形状为(总Z层数, res_y, res_x, 6)。最后一维是[Ex, Ey, Ez, Hx, Hy, Hz]六个场分量。你可以方便地使用matplotlib等工具将场分布可视化这对于理解器件内部的光场增强、模式耦合等物理现象至关重要。4. 深入核心数值实现的关键技术与避坑指南4.1 卷积矩阵生成精度与效率的权衡RCWA中至关重要的一步是将实空间的介电常数分布ε(x,y)转换到傅里叶空间形成卷积矩阵[ε]和[1/ε]。meent在convolution_matrix.py中实现了三种方法对应不同的fft_type。离散傅里叶级数 (DFS,fft_type0)这是最直接的方法对离散的像素值进行FFT得到傅里叶系数。然而对于突变界面如从空气到硅DFS会引入严重的吉布斯振荡导致收敛缓慢甚至错误。meent提供了一个improve_dft选项默认开启其原理是通过上采样复制像素来增加采样频率从而缓解混叠效应。这被称为“增强型DFS”。在我的测试中开启此选项能将与高精度参考解如Reticolo的误差降低约3个数量级。连续傅里叶级数 - 栅格 (CFS-raster,fft_type1)这种方法将每个像素视为一个均匀的矩形块并对其进行解析傅里叶变换。对于每个谐波(m, n)其傅里叶系数ε_{m,n}可以通过二维sinc函数精确计算。这本质上是对阶梯近似结构的精确傅里叶分析精度远高于DFS是meent中精度最高的栅格建模选项。但请注意它目前不支持自动微分。连续傅里叶级数 - 矢量 (CFS-vector,fft_type2)这是精度最高的方法。它直接对解析的几何形状如矩形、椭圆进行傅里叶变换。meent会将旋转后的形状分解为多个与坐标轴对齐的小矩形然后对每个小矩形进行解析积分并求和。n_split参数控制分解的精细程度值越大对原始形状的近似越好计算量也越大。避坑指南傅里叶截断阶数 (Fourier Order) 的选择这是RCWA计算中最重要的参数之一直接决定了计算量和精度。选择过小结果不收敛选择过大计算时间激增复杂度 ~ O(N^3)且可能引入数值误差。经验法则对于折射率对比度不高的结构如硅/空气通常取period/wavelength * 10作为起始点进行扫描。例如周期1000nm波长900nm可以从order10开始。收敛性测试必须进行。针对你的具体结构逐步增加fourier_order观察关心的衍射效率如1级透射是否趋于稳定。当连续增加阶数结果变化小于你的容忍度如1e-3时即可认为收敛。矢量建模的额外考虑在矢量建模中由于边界是光滑的通常比栅格建模收敛得更快可以用更低的阶数获得相同精度。内存警告傅里叶阶数N会生成(2N1)^2大小的矩阵2D光栅。N100时矩阵维度是201x201进行特征值分解和矩阵求逆操作对内存和算力要求很高。务必根据你的硬件条件合理设置。4.2 增强透射矩阵法 (ETM)攻克数值不稳定性在标准的TMM实现中如原文公式61需要计算包含X_ℓ矩阵逆的项。X_ℓ是一个对角矩阵其元素为exp(-k0 * q * d)其中q是特征值。对于厚层或衰减很快的倏逝波模式q的实部很大导致exp(-k0 * q * d)的值极其微小甚至低于机器精度。求逆这样的矩阵会引发严重的数值不稳定结果可能溢出或产生巨大误差。meent采用了增强透射矩阵法ETM来解决这一问题。其核心思想是避免直接对X_ℓ求逆。如原文第G.5节所示通过巧妙的变量代换T A^{-1} X T将问题转化为从最后一层向前一层的递推关系在递推过程中X_ℓ只参与乘法运算而不需要求逆。这极大地提升了数值稳定性使得仿真厚层结构或高折射率对比度结构成为可能。在实际使用中你无需手动操作meent的_solve()方法内部已经实现了ETM。但了解这一原理非常重要因为当你自己编写RCWA代码或使用其他某些库时可能会遇到“矩阵奇异”的报错其根源往往就在这里。4.3 特征值问题的求解后端性能深度剖析RCWA中最耗时的步骤是求解广义特征值问题对于各向同性介质可化为标准特征值问题。meent支持三种后端其性能差异显著选择不当会极大影响工作效率。我基于一个固定的1D光栅结构在三种不同配置的服务器上对64位和32位精度下的NumPy、JAX、PyTorch后端在CPU和GPU上的性能进行了全面测试。以下是我的核心发现和建议1. 后端选择策略追求极致稳定与兼容性选NumPy (CPU)。它的线性代数库如OpenBLAS或MKL经过多年优化非常稳健是交叉验证结果的“金标准”。缺点是慢且不支持AD。进行梯度优化研究首选PyTorch。它的生态完善GPU支持好自动微分接口直观。在我的测试中PyTorch CPU版本在大多数情况下比NumPy快1.5-2倍。PyTorch GPU版本在中小规模矩阵运算上也能带来加速但对于最大的特征值分解Eig这一步由于PyTorch/JAX目前仍调用CPU的LAPACK库数据在CPU/GPU间传输会产生开销导致加速比并不理想有时甚至更慢。尝试JAXJAX的JIT编译理论上能带来性能提升但在处理RCWA这种动态形状矩阵维度随fourier_order变化的问题时每次参数变化都可能触发重新编译产生显著开销。在测试中JAX在处理大阶数FTO200时经常因内存问题失败。目前不建议用于生产环境的大规模计算除非你的问题规模固定且可完全静态编译。2. 精度选择64位 vs 32位32位 (float32/complex64)计算速度通常更快内存占用减半。对于许多光学仿真如果设计目标效率在1%量级32位精度可能足够。64位 (float64/complex128)强烈推荐用于最终验证和发表结果。32位精度只有约7位有效数字在计算大量级数求和、处理非常小或非常大的特征值时累积误差可能不可忽视。特别是当你的优化目标涉及极高效率如99.9%时32位误差可能导致错误结论。一个关键陷阱NumPy的eig函数在输入32位数组时内部计算仍以64位进行最后将结果截断为32位返回。因此NumPy的32位模式并不会节省特征值分解的计算时间节省的只是内存拷贝和后续矩阵运算的时间。而JAX和PyTorch的32位eig是真正的32位运算。3. 硬件利用心得CPU核心数RCWA的矩阵运算特征值分解、矩阵乘法是高度并行的。确保你的NumPy/PyTorch链接了多线程BLAS库如OpenBLAS, MKL并设置合适的线程数如export OMP_NUM_THREADS8。GPU的适用场景当前GPU的加速收益主要体现在卷积矩阵生成、场分布重建等可并行操作上。对于最耗时的特征值分解由于缺乏专用的GPU实现瓶颈仍在CPU。因此不要期望换上GPU就能获得数量级的提升。只有当你的工作流包含大量重复的、小规模矩阵运算如优化迭代中前向计算时GPUPyTorch的组合才可能显示出优势。下表总结了不同场景下的后端选择建议应用场景推荐后端精度关键理由结果验证、交叉比对NumPy (CPU)64位最稳定作为基准参考梯度优化研究PyTorch (CPU/GPU)64位(验证)/32位(迭代)AD支持好生态成熟CPU版本已较快大规模参数扫描无梯度NumPy (CPU)32位/64位稳定内存管理可靠探索性代码开发JAX (CPU)64位语法简洁适合快速原型5. 常见问题排查与实战技巧5.1 仿真结果不收敛或异常这是最常见的问题。请按照以下清单逐步排查检查傅里叶阶数这是首要怀疑对象。逐步增加fourier_order观察结果是否趋向一个稳定值。如果结果剧烈振荡或不收敛很可能阶数太低。检查网格分辨率栅格建模对于fft_type0确保ucell的像素数足够多能够准确描述结构特征。一个经验法则是每个波长范围内至少需要10-20个像素。可以尝试增加像素数或切换到fft_type1(CFS-raster) 看结果是否变化。检查材料参数确认折射率是复数n ik吗k是消光系数代表吸收。如果材料有吸收k0所有衍射级次效率之和应小于1。如果k0但效率和大于1肯定出错了。检查边界条件确认n_I上覆层和n_II衬底设置正确。特别是当光从衬底入射时不要弄反。验证能量守恒计算所有传播级次k_z为实数的反射效率sum(de_ri)和透射效率sum(de_ti)之和。在无吸收的情况下它应该非常接近1例如0.999~1.001。如果偏差很大1%说明计算可能有问题。与简单案例对比仿真一个已知解析解或结果非常简单的结构如均匀介质层应只有0级反射和透射或一个非常浅的光栅可用微扰论验证。用meent计算看结果是否合理。5.2 自动微分AD梯度计算相关问题meent支持通过JAX和PyTorch后端进行自动微分这是其核心优势之一。梯度为None或错误确保参数requires_gradTrue在PyTorch中需要将优化变量用torch.tensor(..., requires_gradTrue)定义。检查计算图确保整个仿真流程从结构参数到最终目标函数如衍射效率都在同一个后端框架内完成。避免中途混用NumPy数组。使用正确的损失函数AD计算的是标量损失函数相对于输入参数的梯度。确保你的目标如-de_ti[1]表示最大化1级透射是一个标量。梯度爆炸或消失这在拓扑优化中很常见。像素的介电常数作为变量梯度可能非常大。尝试使用梯度裁剪torch.nn.utils.clip_grad_norm_或更小的学习率。考虑在优化中引入平滑滤波或投影步骤以稳定优化过程并获得更物理可实现的结构。形状导数 vs 拓扑导数拓扑导数对应于栅格建模 (fft_type0)。梯度是目标函数对每个像素介电常数的导数。优化会自由地改变每个像素的材料拓扑连通性不固定。形状导数对应于矢量建模 (fft_type2)。梯度是目标函数对几何参数如长度、半径、角度的导数。优化过程中结构的拓扑保持不变。根据你的制造约束是灰度曝光还是限定几何形状选择合适的优化方式。5.3 性能优化技巧预热运行对于JAX和PyTorch第一次运行会因为编译或图构建而较慢。在正式计时或优化循环开始前先使用一组典型参数运行几次仿真让系统完成编译。批量计算如果需要扫描波长或角度尽量使用向量化操作。meent的部分函数支持批量输入或者你可以用for循环结合vmap(JAX) 或torch.vmap(PyTorch) 来实现并行计算这比串行循环快得多。控制内存大傅里叶阶数会消耗巨量内存。如果遇到内存错误OOM首先尝试降低fourier_order到刚好收敛的程度。其次可以尝试使用32位精度。在PyTorch中使用torch.cuda.empty_cache()及时清理GPU缓存。选择性计算如果只关心反射或透射meent的conv_solve是最高效的。只有需要分析场分布时才调用calculate_field因为场重建的计算量和内存消耗也很大。6. 案例一维超光栅衍射器优化让我们用一个简化的案例来串联所有知识点。目标是设计一个位于二氧化硅衬底上的一维硅超光栅在正入射TM偏振光波长900nm下将光高效地衍射到1级透射方向。import torch import meent import numpy as np import matplotlib.pyplot as plt # 1. 初始化仿真器 (使用PyTorch后端以支持AD) period 1000 # nm wavelength 900 # nm n_silicon 3.48 0.001j # 硅在900nm附近的复折射率微小吸收 n_sio2 1.45 fourier_order 15 # 初步测试后确定的收敛阶数 # 2. 定义可优化结构 (栅格建模拓扑优化) num_pixels 64 # 初始化设计变量用一个可训练的Tensor表示每个像素的“材料密度”范围[0,1] rho torch.rand(num_pixels, requires_gradTrue) * 0.5 0.25 # 初始值在0.25-0.75之间 rho torch.nn.Parameter(rho) # 将其注册为模型参数 # 3. 定义优化器 optimizer torch.optim.Adam([rho], lr0.05) # 使用Heaviside投影将连续密度二值化并加入平滑滤波避免棋盘格效应 beta 1.0 eta 0.5 def projection(x, beta, eta): 可微分的二值化投影函数 return torch.tanh(beta * eta) torch.tanh(beta * (x - eta)) / torch.tanh(beta * (1 - eta)) # 4. 优化循环 num_iterations 200 efficiency_history [] for i in range(num_iterations): optimizer.zero_grad() # 将密度变量投影到[0,1]并二值化 rho_proj projection(rho, beta, eta) # 可选应用卷积滤波进行平滑 # rho_filtered ... # 根据二值化密度生成介电常数分布硅或空气 epsilon rho_proj * (n_silicon**2 - 1.0) 1.0 # 空气介电常数为1 # 构建单元结构 (单层1D) ucell epsilon[None, None, :] # 形状变为 (1, 1, num_pixels) # 创建仿真实例 mee meent.call_mee( backend2, # PyTorch grating_type0, pol1.0, # TM偏振 (pol1 - ψ0°) n_I1.0, # 上覆空气 n_IIn_sio2, # 衬底二氧化硅 theta0.0, phi0.0, wavelengthwavelength, fourier_orderfourier_order, period[period], ucellucell.detach().cpu().numpy(), # meent接收numpy数组 thickness[300], # 光栅层厚度300nm fft_type0, improve_dftTrue, ) # 运行仿真 de_ri, de_ti mee.conv_solve() # 目标最大化1级透射效率 target_efficiency de_ti[fourier_order 1] # 索引对应1级 loss -target_efficiency # 最小化负效率 # 反向传播 loss.backward() optimizer.step() # 记录并打印 efficiency_history.append(target_efficiency.item()) if i % 20 0: print(fIter {i:3d}, Efficiency: {target_efficiency.item():.4f}, Loss: {loss.item():.4f}) # 逐渐增加二值化强度 if i % 40 0 and i 0: beta min(beta * 1.2, 50) # 缓慢增加beta使投影更接近阶跃 # 5. 结果可视化 plt.figure(figsize(12, 4)) plt.subplot(1, 3, 1) plt.plot(efficiency_history) plt.xlabel(Iteration) plt.ylabel(1st Order Transmission Efficiency) plt.grid(True) plt.subplot(1, 3, 2) final_structure projection(rho, beta, eta).detach().cpu().numpy() plt.bar(range(num_pixels), final_structure) plt.xlabel(Pixel Index) plt.ylabel(Silicon Density (0Air, 1Si)) plt.title(Optimized Binary Structure) # 6. 验证最终性能 (用高精度模式再算一次) mee_verify meent.call_mee( backend0, # 换用更稳定的NumPy验证 grating_type0, pol1.0, n_I1.0, n_IIn_sio2, theta0.0, phi0.0, wavelengthwavelength, fourier_order30, # 使用更高的阶数验证收敛性 period[period], # 将连续密度二值化为0或1 ucell(final_structure 0.5).astype(float)[None, None, :] * (n_silicon**2 - 1) 1, thickness[300], fft_type1, # 使用更高精度的CFS-raster ) de_ri_final, de_ti_final mee_verify.conv_solve() final_eff de_ti_final[30 1] # 1级 plt.subplot(1, 3, 3) orders range(-5, 6) # 显示-5到5级 plt.bar(orders, de_ti_final[30 np.array(orders)], alpha0.7, labelTransmission) plt.bar(orders, de_ri_final[30 np.array(orders)], alpha0.7, labelReflection, bottomde_ti_final[30 np.array(orders)]) plt.xlabel(Diffraction Order) plt.ylabel(Efficiency) plt.legend() plt.title(fFinal Spectrum, 1 order T {final_eff:.3f}) plt.tight_layout() plt.show()这个案例展示了使用meent进行拓扑优化的完整流程定义变量、构建仿真、计算目标、反向传播、更新结构。其中包含了可微分投影、平滑滤波等实用技巧来获得物理上更合理的二值化结构。最后我们用更高精度更高傅里叶阶数和CFS方法验证了优化结果的性能。通过这个从原理到实践从单次仿真到优化设计的完整旅程我们可以看到meent不仅是一个功能强大的RCWA仿真器更是一个连接物理模型、数值计算和优化算法的研究平台。理解其背后的原理能帮助你在遇到问题时快速定位掌握其使用的技巧则能让你在光子器件设计的探索中更加得心应手。