三维各向异性光子晶体的快速仿真算法.pdf

第 卷第 期 年 月东 南 大 学 学 报(自 然 科 学 版)().:./.三维各向异性光子晶体的快速仿真算法杨浩楠 吕星龙 李铁香 杨占山(东南大学数学学院 南京)(南京应用数学中心 南京)(青海民族大学数学与统计学院 西宁)摘要:为获得三维各向异性光子晶体的带隙基于 网格设计了一套求解其能带结构的快速仿真算法.首先 采用有限差分方法对 方程组进行离散 通过结合早期在 氏网格上的工作 对离散所得的 特征值问题的系数矩阵结构进行分析 给出其显式奇异值分解 并利用零空间压缩方法给出无零空间的标准特征值问题形式 在结合求逆 方法和共轭梯度法以及利用快速傅里叶变换大幅加速系数矩阵与向量乘法的基础上 设计出针对三维各向异性光子晶体能带结构的快速数值仿真算法.数值实验表明相比于商业软件 该套算法不仅数值结果准确迭代算法所需的平均次数低于 次且总计算时间少于.展现了算法在结合图形处理单元()高性能计算技术后的有效性与高效性.关键词:光子晶体 网格 各向异性介质 零空间压缩方法 快速傅里叶变换中图分类号:.文献标志码:文章编号:()()()():.().().:()收稿日期:.作者简介:杨浩楠()男硕士生吕星龙(联系人)男博士 .基金项目:国家自然科学基金资助项目().引用本文:杨浩楠吕星龙李铁香等.三维各向异性光子晶体的快速仿真算法.东南大学学报(自然科学版)():.:./.:/.光子晶体是由空间中多种介质材料周期排列而成的人工微结构材料.光子晶体的概念于 年由 和 各自独立提出.光子带隙是光子晶体的基本特性 数值上可通过求解能带结构获得 因此对光子晶体的能带结构设计快速的仿真算法 在理论研究和应用设计上都具有重要意义.光在光子晶体中的传播通常由 方程组来描述 计算光子晶体的能带结构本质上是对应于求解一系列广义特征值问题的数个最小正特征值以及相应的特征向量.近年来 已经有多种计算光子晶体能带结构的数值方法 包括平面波展开法、有限元方法 以及时域有限差分方法 等.平面波展开法通常用于处理周期性边界条件的 方程组 计算大量特征对的时空复杂度较高.有限元方法和有限差分方法可用于处理各种类型的偏微分方程和边界条件 对三维光子晶体来说 通常离散后的特征值问题规模庞大 这对数值算法的高效性提出了巨大挑战.文献利用 氏有限差分方法提出针对各向同性三维光子晶体的一套快速算法.当材料介质为各向异性时 由于介电常数不再为标量函数而表现为矩阵值函数 氏有限差分方法将不再适用.本文将对三维各向异性光子晶体的能带结构设计一套基于 网格的快速数值仿真算法.氏有限差分方法是著名的专门处理电磁场问题的离散方式 特点为电场分量位于网格的边中心 磁场位于面中心 由此形成电场与磁场的 分量环绕形式.网格是由 氏网格及其 个平移网格叠加而成的电场或磁场的 个分量在同一点上均有定义 因此相较于 氏网格更适合用于处理各向异性介质.利用 网格和有限差分方法 对光子晶体的 特征值问题进行数值离散 即可得到相应的广义特征值问题.通过对离散特征值问题系数矩阵的结构分析 提出一种针对计算最小正特征值及相应特征向量的数值算法.本文算法适用于三维光子晶体的所有 种布拉维晶格结构 并结合了快速傅里叶变换()和图形处理单元()等技术是一种有效且高效的数值仿真算法.数学模型三维光子晶体通常由如下频域无源的 方程组来描述:()()()()()()()()()()式中 表示频率()、()、()、()分别表示位于 ()处的电场强度、磁场强度、电位移矢量和磁感应强度.同时 对于线性非色散介质 和 满足如下本构关系:()()()()()()()()式中()、()分别为材料的介电常数和磁导率.由光子晶体的周期性结构特点可知介电常数()和磁导率()满足如下周期性条件:()()()()式中 为晶格平移向量.本文考虑各向异性介电质材料 磁导率()介电常数()为 正定矩阵.对于该张量形式的介电常数()本构关系()处的乘积()()要求()的 个方向分量的离散网格点定义在同一网格点处如图()所示.而在标准的 氏网格上 个方向分量定义在不同的网格点处如图()所示不适用于各向异性材料.故本文将采用由 氏网格(子网格)及其 个平移网格(子网格 )组合而成的 网格对式()和()进行数值离散.利用式()和()消去()便可以得到如下 特征值问题:()()()()()()()式中 为待求的特征值.由于光子晶体的周期性结构特点 该特征值问题只需在单个晶胞上进行求解 且根据布洛赫定理 电场()满足如下拟周期性条件:()()()式中 为第一布里渊区的布洛赫波向量.该性质同样适用于磁场()、电位移矢量()以及磁感应强度().在利用 网格对 特征值问题()进行离散时 散度为零条件可自动满足 因此 特征值问题可离散为如下广义特征值问题:东南大学学报(自然科学版)第 卷:/.()网格 ()子网格 ()子网格()子网格 ()子网格 图 网格以及 和 的分量配置 式中 为双旋度算子 的离散矩阵 为介电常数()的离散矩阵 为离散()后向量化的结果.在本文中表示 维的单位矩阵 和 分别表示矩阵的转置和共轭转置 为 积 ()表示由矩阵、构成的分块对角阵()表示将矩阵向量化处理.网格上的数值离散.和 的离散在利用有限差分方法进行离散时 需预先对三维光子晶体的原胞(见图 中的平行六面体区域)进行裁剪和粘贴 从而拼接成计算晶胞(见图 中的长方体区域)对所有 种布拉维晶格的详细裁剪方式请参考文献.记 沿、三个方向上的长度为、.图 三维光子晶体原胞和计算晶胞示意图对 在、三个方向上分别均匀剖分、份 并记、为相应的步长()表示 网格上的任意网格点()简化表示为().而中点记为().如图 所示 网格由 氏网格(子网格)及其 个平移网格(子网格 )组合而成平 移 量 分 别 为、()及 .因此其上的()与()的分量按以下 部分组成.首先在图()的子网格 上记()和()的分量分别为()、()、()和()、()、().类似地 图()、()、()上()与()的分 量 则 可 以 通 过 平 移、得到.将()与()在、三个方向上的分量在上述 个网格上的取值分别记为、以及、.以子网格 为例定义第 期杨浩楠等:三维各向异性光子晶体的快速仿真算法:/.(:)(:)(:)(:)(:)(:)子网格、上的离散电磁场、与、都可类似以上形式定义.最 后 用 、来表示如图()所示的电磁场离散在整个 网格上电磁场的离散结果.介电常数 及磁导率 的离散对于三维各向异性光子晶体定义()(:)()(:)()(:)()(:)式中 ()、()、()和()分别为 的分量 在 网格中沿 方向的边中心点、沿 方向的边中心点、沿 方向的边中心点和体中心点上的取值.构造 矩阵 有()()()()()()()()()()()()()()()()()()()()()()()()此 矩阵即为介电常数 在 网格上的数值离散.磁导率 的离散矩阵可类似得到.由于本文考虑()因此 的离散矩阵为.定理 矩阵 是正定矩阵.证明:定义 其中 利用 对矩阵 进行重排 得到()式中()()()().这里、和 分别对应介电常数 在 网格中沿、方向以及在体中心点上的数值离散形式.由介电常数 的正定性可知、和 均为正定矩阵 因此矩阵 亦正定.证毕.网格上 的离散在本节中 将推导 网格上式()中的双旋度算子的显式矩阵表示.先定义如下的矩阵函数:()以体心立方晶格为例子 在 氏网格中 的离散矩阵表示为()其中、的形式如下所示:()()()()()()/()东南大学学报(自然科学版)第 卷:/.()()()/().的矩阵表示为构造 网格下式()的离散形式需要分别对图 中子网格、进行考虑.以子网格 为例分别考虑 在 氏网格中面心点()、()以及()上的离散.由式()并结合本构关系式()有()()()()()()()()()()()()()()()因此子网格 上式()的离散矩阵表示为()()同理 子网格、上式()的离散矩阵分别表示为()()()()()()综上 式()在 网格上的离散矩阵表示为 ()其中 ()()()()().的矩阵表示类似.节中的讨论 结合本构关系式()以及.节 式()在图 中 个子网格上的矩阵分别表示为()()()()()()()()()()()()因此 网格上式()的矩阵表示为 ()由于 网格本身满足散度为零条件()式()的离散格式可以通过联立式()和()得到即 ()其中.至此 三维各向异性光子晶体能带结构的计算即转化为求解一系列以波向量 为自变量的广义特征值问题()的几个最小正特征值.快速仿真算法基于前期对各向同性三维光子晶体的研究对 进行奇异值分解并在此基础上设计求解式()中特征值问题的快速算法.的奇异值分解考虑到 是块对角矩阵 由()、()、()、()四部分组成只需对这 部分分别进行奇异值分解 再按照式()中的形式组装即可得到 的奇异值分解.由文献可知存在酉矩阵 可将式()中定义的、同时对角化 即 式中为 的特征值构成的对角矩阵.具有 结构 它们的具体形式参见文献.因此 矩阵()、()、()、()均存在如下形式的矩阵分解:()()其中 ().()分别取值()、()、()及().基于此 定义 ()()()第 期杨浩楠等:三维各向异性光子晶体的快速仿真算法:/.式中()().则有 的奇异值分解(/)为了 清 楚 起 见 分 别 记 种 不 同 情 况 的为 以方便引出如下 的奇异值分解.定理(的奇异值分解)存在酉矩阵()()以 及 它 们 的 前 列 构 成 的 矩 阵()()使得()(/)(/)()(/)(/)()(/)(/)()(/)(/)令 ()()(/)则有 的奇异值分解 ()根据式()和 即可有 ()由式()可知离散矩阵 具有/维度的零空间 这将会严重减慢计算最小正特征值时迭代算法的收敛速度.接下来将利用定理 和零空间压缩方法 去除特征值问题()中的零特征值.无零空间的标准特征值问题令 .结合式()重写式()可将该广义特征值问题转化为无零空间的标准特征值问题即 ()式中.由此对三维光子晶体能带结构的计算转化为求解式()中的数个最小正特征值.对式()的数值计算有以下好处:问题的维度由原先的 降为了 且保留了所有的非零特征值 更重要的是 没有了零特征值的干扰 求逆 方法可轻松用来求解数个最小正特征值.在对式()利用求逆 方法计算的过程中 需要求解形如下式的线性系统:()在实际计算中 将.节中 矩阵()替换为它的逆()采用与矩阵 相同的构造方式即可显式构造出矩阵 大大减少了求解式()的计算量.又有 为正定矩阵且其条件数 ()较小()()因此可以直接利用共轭梯度法(方法)进行求解.利用 方法求解式()的主要计算量在于矩阵 与向量的乘法.根据定理 中矩阵 的定义 与向量的乘法运算实际为对角矩阵与向量的乘法和稠密阵 与向量的乘法.而矩阵 与向量的乘法可以利用 代替既避免了稠密阵 的存储 也极大地加速了矩阵向量乘法.综上 对广义特征值问题()中系数矩阵的奇异值分解 结合压缩零空间技术以及基于 的快速矩阵向量乘法 使得可以对三维各向异性光子晶体能带结构计算问题设计如下快速算法.算法 三维各向异性光子晶体的快速算法输入:晶体位移向量、材料的介电常数(包括内部介质介电常数 及背景材料介电常数)网格数、和波向量 所需特征值的个数.输出:式()的 个目标特征对()(、分别表示第 个特征值与特征向量).对、根据.节中的讨论构造矩阵 对、按照定理 中格式生成、利 用 求 逆 方 法 求 解 式()特 征 对().在每次迭代中 利用矩阵向量乘法 和 结合定理 中 的定义利用 方法求解线性系统()计算 计算 数值实验在本节中针对算法 开展相应的数值实验.由于算法 中主要运算均为向量运算以及傅里叶东南大学学报(自然科学版)第 卷:/.变换运算 因此均可利用 高性能计算技术进行大幅加速.本节的数值实验均在配置有 的服务器上完成.考虑体心立方晶格结构的三维光子晶体 其平行六面体的原胞结构如图()所示 拼接后形成的长方体计算晶胞结构如图()所示 具体拼接细节参考文献.令晶格常数 .体心立方晶格的晶格平移向量分别为 .在该模型中 介质形状均为圆柱且圆柱的半径 .背景介质为空气.磁导率为 ()介 电 质 的 介 电 常 数 .对图()中计算晶胞沿、三个方向的网格剖分数分别为、.波向量 的个数为.此时式()中的广义特征值问题的矩阵维度为 式()中的无零空间的标准特征值问题的矩阵维度为 .()体心立方晶格原胞()体心立方晶格计算晶胞图 体心立方晶格示意图在计算光子晶体的能带结构时 对于每一个波向量 计算无零空间标准特征值问题()的 个最小正特征值 取 子空间的维度为.需要说明的是 网格由 个子网格耦合而成 因此离散特征值问题会存在 倍的增根 即 条能带变为 条相近的能带 可取其均值作为计算的结果.为了清楚展示本文所提仿真算法的数值效果将 个特征值的结果与商业软件 所计算的能带结构(个特征值)进行对比 如图 所示.图中归一化频率为/()为真空中的光速.从图中可以看出本文数值结果与 的计算结果十分吻合 因此算法是可靠的.图 能带结构计算结果本文所提算法可以充分发挥 高性能计算技术的优势每一个波向量 所对应的计算时间以及算法 中求逆 方法和 方法的迭代次数如图 所示.图()为 与本文算法在能带计算上所花费的时间其中 软件中离散网格单元数约为 (本文的有限()计算时间()求逆 方法迭代次数()方法的平均迭代次数图 计算时间及迭代算法次数第 期杨浩楠等:三维各向异性光子晶体的快速仿真算法:/.差分网格数约为 ).从中可以看到对于每一个波向量的计算 平均需要花费 而本文算法均不到 完成整个能带计算需要大约 而本文算法在 加速下不超过.从图()可知对一个波向量求解 个最小正特征值 求逆的 方法的总迭代次数不超过.图()则说明平均每步 方法迭代中 利用 求解线性系统()的迭代次数均不超过.该实验说明基于 高性能计算技术算法 对于求解式()中的大规模矩阵特征值问题是十分高效的.结论)本文利用 网格和有限差分方法对三维各向异性光子晶体的 特征值问题进行数值离散.)在充分挖掘离散系数矩阵结构的基础上给出显式的奇异值分解公式 并将离散广义特征值问题压缩为无零空间的标准特征值问题.该过程去除了/的零特征值 有效地实现了离散模型的降阶 同时系数矩阵的显式奇异值分解所蕴含的 结构使得快速求解大规模特征值问题成为可能.)结合求逆 方法和 方法设计了针对最小正特征值以及相应特征向量的算法 即完成了求解三维光子晶体能带结构的快速仿真算法.)基于迭代算法中良态的线性系统和系数矩阵向量乘法中蕴含的 结构结合 高性能计算技术实现了该算法 相关的仿真软件包发布在 (.).)数值实验部分针对同一算例与商业软件 相比本文算法不但数值结果与之吻合且 及 两个迭代算法的平均次数分别低于 与 总计算时间少于.显示了算法在结合 后的有效性与高效性.参考文献().():.:./.():.:./.():.:./.():.:./.():.:./.():.:./.:.:./.():.:./.:.:./.():.:./.():.:./.():.:./.:.:./.:.():.:./.东南大学学报(自然科学版)第 卷