CN107247686B - 一种基于并行算法的fetd仿真模拟方法 - Google Patents
一种基于并行算法的fetd仿真模拟方法 Download PDFInfo
- Publication number
- CN107247686B CN107247686B CN201710361396.6A CN201710361396A CN107247686B CN 107247686 B CN107247686 B CN 107247686B CN 201710361396 A CN201710361396 A CN 201710361396A CN 107247686 B CN107247686 B CN 107247686B
- Authority
- CN
- China
- Prior art keywords
- equation
- time
- fetd
- solution
- time advance
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
- G06F17/12—Simultaneous equations, e.g. systems of linear equations
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F30/00—Computer-aided design [CAD]
- G06F30/20—Design optimisation, verification or simulation
- G06F30/23—Design optimisation, verification or simulation using finite element methods [FEM] or finite difference methods [FDM]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F2111/00—Details relating to CAD techniques
- G06F2111/10—Numerical modelling
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Mathematical Physics (AREA)
- Theoretical Computer Science (AREA)
- Computational Mathematics (AREA)
- Mathematical Analysis (AREA)
- Mathematical Optimization (AREA)
- Pure & Applied Mathematics (AREA)
- General Engineering & Computer Science (AREA)
- Data Mining & Analysis (AREA)
- Operations Research (AREA)
- Software Systems (AREA)
- Databases & Information Systems (AREA)
- Algebra (AREA)
- Computer Hardware Design (AREA)
- Evolutionary Computation (AREA)
- Geometry (AREA)
- Complex Calculations (AREA)
Abstract
本发明公开了一种基于并行算法的FETD仿真模拟方法,属于计算电磁学领域。该方法在现有算法的基础上,通过并行计算的方式替代常规的串行计算,来实现对方程循环迭代的求解。在实现并行求解的过程中,通过引入矩阵对方程组变换来克服时间推进方程中三个相邻时刻未知量之间的依赖关系,从而为实现并行提供可行性。该方法在保证计算精度的基础上能够解决FETD中的时间推进方程循环迭代求解效率低下的问题。
Description
技术领域
本发明属于计算电磁学领域,具体涉及一种基于并行算法的FETD仿真模拟方法。
背景技术
自从1864年Maxwell方程组被提出,人们对电磁波的研究不断深入,电磁理论的应用已经遍及生命科学、医学、材料科学和信息科学等其他科学领域。对于复杂的电磁问题,通过实验方法往往会受到试验测试成本过高研究周期长、甚至试验无法实现等问题的限制;而理论分析时由于理论模型十分复杂,经典电磁学解析求解时存在较大的局限性甚至根本无法实现。为了解决这类电磁问题,随着计算机软硬件计算的发展,结合电磁理论和计算数学各种数值计算方法相继被提出,计算电磁学这门交叉学科也应运而生,计算电磁学已成为现代计算电磁理论不可或缺的一部分。
时域有限元方法(FETD)作为计算电磁学领域中的一种时域数值计算方法,它既能对复杂几何结构进行模拟计算,又能通过对计算结果进行离散傅里叶变换而得到结构的宽频带特性。FETD继承了频域有限元方法的有点的同时,它还能直接在时域内进行计算,因而FETD在计算电磁学领域得到深入的发展和广泛的应用。FETD从麦克斯韦方程组出发通过插值基函数将未知量展开来进行空间离散,并且通过稳定的时间差分格式来进行时间离散,从而实现对电磁问题的数值求解。在运用FETD将待求问题的微分控制方程进行空间离散和时间离散后可以得到待求问题的时间推进方程,通过时间推进方程可以从初值时刻未知量的值推导出后面任意时刻未知量的值。对于运用FETD求解电磁问题的过程,[The FiniteElement Method in Electromagnetics,529-577页,作者:J.M.Jin]一文中有详细的介绍。这个过程需要不断的循环迭代求解,随着循环次数的增多,该循环迭代求解过程在编程实现上需要消耗大量的时间从而影响计算效率。为了解决这个难题,设计了一种并行优化计算的方法来避免对时间推进方程的循环迭代求解,以达到提高计算效率的目的。
发明内容
本发明的目的是提供一种方法来解决FETD计算过程中的时间推进方程循环迭代求解效率低下的问题。该方法通过并行处理能有效避免循环求解过程,从而提高计算效率。
为了实现上述目的,本发明的技术方案是:一种基于并行算法的FETD仿真模拟方法,包括以下步骤:
A.确定需要分析的电真空器件结构;
B.对步骤A的器件结构进行建模,建立对应的几何结构模型;
C.确定电真空器件结构的电磁学边值问题的控制微分方程形式;
D.采用四面体单元网格剖分求解区域;
E.用插值基函数将控制微分方程中的待求未知量进行空间离散展开,并运用标准变分原理得到边值问题关于时间偏微分的有限元方程组;
F.选择稳定的时间差分格式(如中心差分格式、newmark-β差分格式)对步骤E中的有限元方程组进行时间离散,得到边值问题的时间推进方程。
G.采用并行算法计算步骤F中的时间推进方程的迭代求解过程。
与现有技术相比,本发明的有益效果:利用本发明提出的一种基于FETD的时间推进方程迭代求解过程的并行实现算法,在保证计算精度的同时能够有效提高该过程计算效率。
附图说明
图1是矩形波导结构网格离散后的示意图。
图2是本发明基于并行算法的FETD仿真模拟方法的流程图。
图3是FETD的时间推进方程迭代求解过程串行实现流程图。
图4是FETD的时间推进方程迭代求解过程并行实现流程图。
具体实施方式
下面结合附图和具体实施方式来详细描述本发明的技术方案。
参照流程图2,一种基于并行算法的FETD仿真模拟方法包括以下步骤:
A.确定需要分析的电真空器件结构。
选取需要分析的电真空器件结构,如波导、矩形窗、盒形窗、同轴窗等,本实施例选取矩形波导。
B.对步骤A的器件结构进行建模,建立对应的几何结构模型;
本发明基于矩形波导结构而言,网格离散后的几何结构模型如图1所示。
C.确定电真空器件结构的电磁学边值问题的控制微分方程形式
对于矩形波导结构而言,对该结构电磁传输特性的分析是在一定的空间内求解Maxwell方程组,该问题的求解域为波导壁界定的空间。从Maxwell方程组出发,在无外加电流激励的情况下可以得到电场E满足的微分方程为:
其中,μ为磁导率,ε为介电常数,σ为电导率,t为时间。为了使公式(1)中的方程解唯一,需要确定波导入射端口面和出射端口面上的边界条件。一般边界条件满足:
D.采用四面体单元网格剖分求解区域。
采用四面体网格剖分求解域是时域有限元方法中的一种公知过程,因此本步骤不再详细描述。剖分后的求解域被人为分割为多个三维四面体网格,从而将连续的几何结构空间转化为离散的网格空间。
E.用插值基函数将控制微分方程中的待求未知量进行空间离散展开,并运用标准变分原理得到边值问题关于时间偏微分的有限元方程组。
如公式(1)所示,对于步骤C中的控制微分方程的求解需要先用插值基函数将待求未知量进行空间离散展开。这个也是公知的过程,因此本步骤不再详细描述,直接给出通过标准变分原理将控制微分方程进行空间离散后的具体形式:
其中{u}为电场按基函数展开的展开系数。各系数矩阵[T]、[Tσ]、[S]、{f}中的矩阵元素满足:
公式(4)中,V是体积,S是边界面面积,Ni和Nj均为插值基函数,插值基函数为矢量棱边基函数时i的取值为:i=1,2,3,4,5,6。同样j的取值范围和i一致。
F.选择稳定的时间差分格式对步骤E中的有限元方程组进行时间离散,得到边值问题的时间推进方程。
在步骤E的基础上,对公式(3)中的方程进行时间离散。这里选取中心差分格式进行离散得到时间推进方程为:
其中,Δt为时间步长,{f}n为激励,{u}n为第n时刻的待求未知量。
G.通过并行算法实现步骤F中的时间推进方程的迭代求解过程。
对于公式(5)中的时间推进方程,在给出激励形式、初始值{u}0和{u}1后,可以求解出{u}2的值;再通过{u}1和{u}2求出{u}3的值,依次类推可以求出后续任意时刻的{u}值。这种方式可以看成一个迭代递推的过程,目前处理该过程采用的方法都是用循环的方式串行实现,其流程如图3所示。流程图3展示了待求问题循环迭代求解的串行过程,初始零时刻时,{u}0,{u}1,f0已知,在满足循环条件的情况下进入循环体可以计算出{u}2的值,接着返回循环条件继续判断是否满足该条件,如果满足循环条件,则继续进入循环体计算{u}3,{u}4,…,{u}n。循环条件是通过迭代次数进行约束,在循环之前需要知道该约束条件。上述过程需要通过不断的循环迭代求解从而得到后续时刻的值,这种串行计算的思想在计算效率上并不能达到最优。
本发明中步骤A-F为常规计算过程,核心发明点为采用并行算法替代现有的串行算法,具体实施过程为:
由公式(5),FETD的时间推进方程可以简化为:
[A]{u}n+1=[B]{u}n+[C]{u}n-1+{f}n (6)
其中[A],[B],[C]为系数矩阵,并且由公式(5)中的各系数矩阵得到。这里对公式(6)两边同乘[A]-1得到:
{u}n+1=[A]-1[B]{u}n+[A]-1[C]{u}n-1+[A]-1{f}n (7)
在主模激励的情况下,{f}n满足:
联立公式(7)和公式(8)可得:
将公式(9)进行简化得到:
{U}n+1=[M]n{U}n (10)
其中:
根据公式(10)并令[M]1=[M]n,[M]2=[M]n[M]n-1,…,[M]n=[M]n[M]n-1…[M]1,则有:
{U}n=[M]1{U}n-1=[M]2{U}n-2=…=[M]n-1{U}1 (11)
也即对公式(3)的求解转化为公式(11)的解。对于公式(11)的迭代求解过程,采用并行算法实现,其流程图如图4所示。流程图4展示了待求问题的并行实现过程,在给定初值的情况下,程序实现时需要创建多个线程,线程的个数根据需要迭代的次数确定。在计算过程中,让多个线程同时工作,例如:线程1计算{U}2的值,线程2计算{U}3的值,线程3计算{u}4的值等等。从流程图4可以看到,各个线程的工作相互独立且同时进行。从加快计算速度和提高计算效率的角度考虑,上述这种并行的处理方式能够解决目前传统方法中存在的问题。这种并行处理迭代求解过程的方式和串行循环处理相比在计算效率上得到优化,能够充分利用计算机的资源。
还需要注意的是对于N×N维矩阵[M],当N取20000时用数组的方式存储该矩阵需要8×20000×20000个字节(3GB)的内存。这种直接存储矩阵元素的方式对内存的消耗过大,在数值计算编程过程中是要避免的。由于矩阵[M]对称正定且矩阵中有很多零元素,而这些零元素是不直接参与计算的,对计算产生作用的只有那些非零元素。从考虑节约内存的角度出发,我们采用稀疏存储的方式来存储矩阵中的非零元素。稀疏存储的方式主要有按行稀疏存储和按列稀疏存储。这里采用按列稀疏存储的方式来进行,对于N阶m个非零元素的矩阵来说,用三个一维数组*ptr,*ind和*val来稀疏存储,ind,val的长度是m,按行顺序记录每个非零元的列标和数值;ptr的长度是N+1,ptr[i]记录第i行第一个非零元的位置。最后一个元素ptr[N]=m。例如对于三维矩阵:
则:
ptr[4]={0,2,4,6};
ind[6]={0,2,1,2,0,2};
val[6]={1.0,2.0,3.0,4.0,5.0,6.0};
因为[M]是按列稀疏方式存储,在计算公式(11)时,可以利用英特尔数学核心函数库(Intel MKL)里的相关库函数进行求解计算。用稀疏的方式存储矩阵并结合Intel MKL里的库函数对公式(11)中的方程进行求解,不仅节约了大量内存空间,更提高了求解效率。
Claims (1)
1.一种基于并行算法的FETD仿真模拟方法,包括以下步骤:
A.确定需要分析的电真空器件结构;
B.对步骤A的器件结构进行建模,建立对应的几何结构模型;
C.确定电真空器件结构的电磁学边值问题的控制微分方程形式;
D.采用四面体单元网格剖分求解区域;
E.用插值基函数将控制微分方程中的待求未知量进行空间离散展开,并运用标准变分原理得到边值问题关于时间偏微分的有限元方程组;
F.选择稳定的时间差分格式对步骤E中的有限元方程组进行时间离散,得到边值问题的时间推进方程;所述时间差分格式为中心差分格式;
G.采用并行算法计算步骤F中的时间推进方程的迭代求解过程;
具体地,所述步骤F中边值问题的时间推进方程为:
其中Δt为时间步长,{f}n为激励,{u}n为第n时刻电场按基函数展开的展开系数,{u}n+1,{u}n-1分别为{u}n相邻时刻的值,各系数矩阵[T]、[Tσ]、[S]、{f}中的矩阵元素满足:
公式(2)中,μ为磁导率,ε为介电常数,σ为电导率,γ为边界等效波导纳,为已知量,V是体积,S是边界面面积;Ni和Nj均为插值基函数,插值基函数为矢量棱边基函数时i的取值为:i=1,2,3,4,5,6,j的取值范围和i一致;
所述步骤G中的并行算法的具体过程为:
由公式(1),FETD的时间推进方程可以简化为:
[A]{u}n+1=[B]{u}n+[C]{u}n-1+{f}n (3)
其中[A],[B],[C]为系数矩阵,由时间推进方程中的各系数矩阵得到,对公式(3)两边同乘[A]-1得到:
{u}n+1=[A]-1[B]{u}n+[A]-1[C]{u}n-1+[A]-1{f}n (4)
在主模激励的情况下,{f}n满足:
联立公式(4)和公式(5)可得:
将公式(6)进行简化得到:
{U}n+1=[M]n{U}n (7)
其中:
根据公式(7)并令[M]1=[M]n,[M]2=[M]n[M]n-1,…,[M]n=[M]n[M]n-1…[M]1,则有:
{U}n=[M]1{U}n-1=[M]2{U}n-2=…=[M]n-1{U}1 (8)
也即对公式(1)的求解转化为公式(8)的解。
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201710361396.6A CN107247686B (zh) | 2017-05-22 | 2017-05-22 | 一种基于并行算法的fetd仿真模拟方法 |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN201710361396.6A CN107247686B (zh) | 2017-05-22 | 2017-05-22 | 一种基于并行算法的fetd仿真模拟方法 |
Publications (2)
Publication Number | Publication Date |
---|---|
CN107247686A CN107247686A (zh) | 2017-10-13 |
CN107247686B true CN107247686B (zh) | 2020-11-17 |
Family
ID=60017417
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN201710361396.6A Active CN107247686B (zh) | 2017-05-22 | 2017-05-22 | 一种基于并行算法的fetd仿真模拟方法 |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN107247686B (zh) |
Families Citing this family (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN107944141A (zh) * | 2017-11-24 | 2018-04-20 | 电子科技大学 | 基于杂交时域间断伽辽金法的时域计算电磁学数值方法 |
CN110909511B (zh) * | 2019-11-22 | 2022-10-14 | 电子科技大学 | 一种无曲面体积分的无粘低速绕流数值模拟方法 |
CN111209675B (zh) * | 2020-01-10 | 2023-04-28 | 南方电网科学研究院有限责任公司 | 电力电子器件的仿真方法、装置、终端设备及存储介质 |
CN113722966B (zh) * | 2021-09-07 | 2022-08-19 | 芯和半导体科技(上海)有限公司 | 一种集成电路板仿真多级分布式并行计算方法 |
Citations (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN101788819A (zh) * | 2010-03-08 | 2010-07-28 | 清华大学 | 大规模生产过程一种基于迭代式分解和流松驰的调度方法 |
CN102609575A (zh) * | 2012-01-19 | 2012-07-25 | 浙江大学 | 一种基于隐式数值积分的电力***暂态稳定仿真方法 |
CN104504174A (zh) * | 2014-11-28 | 2015-04-08 | 南京航空航天大学 | 一种适于伞降空投过程流固耦合数值模拟的自适应网格控制方法 |
CN104578054A (zh) * | 2015-01-12 | 2015-04-29 | 浙江大学 | 一种基于多稀疏矢量道路集的电力***暂态稳定仿真方法 |
CN105911584A (zh) * | 2015-09-25 | 2016-08-31 | 中国科学院地质与地球物理研究所 | 一种隐式交错网格有限差分弹性波数值模拟方法及装置 |
-
2017
- 2017-05-22 CN CN201710361396.6A patent/CN107247686B/zh active Active
Patent Citations (5)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN101788819A (zh) * | 2010-03-08 | 2010-07-28 | 清华大学 | 大规模生产过程一种基于迭代式分解和流松驰的调度方法 |
CN102609575A (zh) * | 2012-01-19 | 2012-07-25 | 浙江大学 | 一种基于隐式数值积分的电力***暂态稳定仿真方法 |
CN104504174A (zh) * | 2014-11-28 | 2015-04-08 | 南京航空航天大学 | 一种适于伞降空投过程流固耦合数值模拟的自适应网格控制方法 |
CN104578054A (zh) * | 2015-01-12 | 2015-04-29 | 浙江大学 | 一种基于多稀疏矢量道路集的电力***暂态稳定仿真方法 |
CN105911584A (zh) * | 2015-09-25 | 2016-08-31 | 中国科学院地质与地球物理研究所 | 一种隐式交错网格有限差分弹性波数值模拟方法及装置 |
Non-Patent Citations (2)
Title |
---|
"时域有限元方法在仿真微波无源器件中的应用";彭聪;《中国优秀硕士学位论文全文数据库 工程科技Ⅱ辑》;20160315(第3期);第C028-14页 * |
"时域有限元电磁计算方法的研究";杜磊;《中国博士学位论文全文数据库 基础科学辑》;20160515(第5期);第A002-18页 * |
Also Published As
Publication number | Publication date |
---|---|
CN107247686A (zh) | 2017-10-13 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN107247686B (zh) | 一种基于并行算法的fetd仿真模拟方法 | |
EP3179415B1 (en) | Systems and methods for a multi-core optimized recurrent neural network | |
WO2010039325A1 (en) | Method for solving reservoir simulation matrix equation using parallel multi-level incomplete factorizations | |
CN104951580B (zh) | 无条件稳定和有条件稳定混合的时域谱元电磁分析方法 | |
Kamranian et al. | An adaptive meshless local Petrov–Galerkin method based on a posteriori error estimation for the boundary layer problems | |
CN109726441B (zh) | 体和面混合gpu并行的计算电磁学dgtd方法 | |
CN113158492B (zh) | 一种时变电磁场的全隐式双时间步计算方法 | |
Haghi et al. | FP-AMG: FPGA-based acceleration framework for algebraic multigrid solvers | |
CN116258042B (zh) | 一种基于ddm的大规模传热异构并行仿真方法 | |
Nakajima | Optimization of serial and parallel communications for parallel geometric multigrid method | |
Dziekonski et al. | GPU-accelerated LOBPCG method with inexact null-space filtering for solving generalized eigenvalue problems in computational electromagnetics analysis with higher-order FEM | |
CN114757347A (zh) | 低比特量化神经网络加速器实现方法和*** | |
CN116150553A (zh) | 一种面向cpu+dcu异构混合架构的稀疏化amg优化方法 | |
CN113239591B (zh) | 面向dcu集群的大规模有限元网格并行分区的方法及装置 | |
Wu et al. | High-performance architecture for the conjugate gradient solver on FPGAs | |
CN111339688B (zh) | 基于大数据并行算法求解火箭仿真模型时域方程的方法 | |
Li | Optimal L2 error estimates for the interior penalty DG method for Maxwell’s equations in cold plasma | |
Carpentieri | Fast iterative solution methods in electromagnetic scattering | |
Kumar et al. | Multi-threaded nested filtering factorization preconditioner | |
Li et al. | ARCHITECT: Arbitrary-precision constant-hardware iterative compute | |
Kansa | Radial basis functions: achievements and challenges | |
Mukaddes et al. | Treatment of Block-Based Sparse Matrices in Domain Decomposition Method | |
Mukaddes et al. | Performance evaluation of domain decomposition method with sparse matrix storage schemes in modern supercomputer | |
Zhang et al. | Accelerating the finite element method using FPGA for electromagnetic field computation | |
Lin et al. | Performance of a massively parallel higher-order method of moments code using thousands of CPUs |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
PB01 | Publication | ||
PB01 | Publication | ||
SE01 | Entry into force of request for substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant |