本發(fā)明屬于計(jì)算電磁學(xué)領(lǐng)域，具體涉及一種基于并行算法的fetd仿真模擬方法。

背景技術(shù)：

自從1864年maxwell方程組被提出，人們對(duì)電磁波的研究不斷深入，電磁理論的應(yīng)用已經(jīng)遍及生命科學(xué)、醫(yī)學(xué)、材料科學(xué)和信息科學(xué)等其他科學(xué)領(lǐng)域。對(duì)于復(fù)雜的電磁問(wèn)題，通過(guò)實(shí)驗(yàn)方法往往會(huì)受到試驗(yàn)測(cè)試成本過(guò)高研究周期長(zhǎng)、甚至試驗(yàn)無(wú)法實(shí)現(xiàn)等問(wèn)題的限制；而理論分析時(shí)由于理論模型十分復(fù)雜，經(jīng)典電磁學(xué)解析求解時(shí)存在較大的局限性甚至根本無(wú)法實(shí)現(xiàn)。為了解決這類電磁問(wèn)題，隨著計(jì)算機(jī)軟硬件計(jì)算的發(fā)展，結(jié)合電磁理論和計(jì)算數(shù)學(xué)各種數(shù)值計(jì)算方法相繼被提出，計(jì)算電磁學(xué)這門交叉學(xué)科也應(yīng)運(yùn)而生，計(jì)算電磁學(xué)已成為現(xiàn)代計(jì)算電磁理論不可或缺的一部分。

時(shí)域有限元方法(fetd)作為計(jì)算電磁學(xué)領(lǐng)域中的一種時(shí)域數(shù)值計(jì)算方法，它既能對(duì)復(fù)雜幾何結(jié)構(gòu)進(jìn)行模擬計(jì)算，又能通過(guò)對(duì)計(jì)算結(jié)果進(jìn)行離散傅里葉變換而得到結(jié)構(gòu)的寬頻帶特性。fetd繼承了頻域有限元方法的有點(diǎn)的同時(shí)，它還能直接在時(shí)域內(nèi)進(jìn)行計(jì)算，因而fetd在計(jì)算電磁學(xué)領(lǐng)域得到深入的發(fā)展和廣泛的應(yīng)用。fetd從麥克斯韋方程組出發(fā)通過(guò)插值基函數(shù)將未知量展開來(lái)進(jìn)行空間離散，并且通過(guò)穩(wěn)定的時(shí)間差分格式來(lái)進(jìn)行時(shí)間離散，從而實(shí)現(xiàn)對(duì)電磁問(wèn)題的數(shù)值求解。在運(yùn)用fetd將待求問(wèn)題的微分控制方程進(jìn)行空間離散和時(shí)間離散后可以得到待求問(wèn)題的時(shí)間推進(jìn)方程，通過(guò)時(shí)間推進(jìn)方程可以從初值時(shí)刻未知量的值推導(dǎo)出后面任意時(shí)刻未知量的值。對(duì)于運(yùn)用fetd求解電磁問(wèn)題的過(guò)程，[thefiniteelementmethodinelectromagnetics，529-577頁(yè)，作者：j.m.jin]一文中有詳細(xì)的介紹。這個(gè)過(guò)程需要不斷的循環(huán)迭代求解，隨著循環(huán)次數(shù)的增多，該循環(huán)迭代求解過(guò)程在編程實(shí)現(xiàn)上需要消耗大量的時(shí)間從而影響計(jì)算效率。為了解決這個(gè)難題，設(shè)計(jì)了一種并行優(yōu)化計(jì)算的方法來(lái)避免對(duì)時(shí)間推進(jìn)方程的循環(huán)迭代求解，以達(dá)到提高計(jì)算效率的目的。

技術(shù)實(shí)現(xiàn)要素：

本發(fā)明的目的是提供一種方法來(lái)解決fetd計(jì)算過(guò)程中的時(shí)間推進(jìn)方程循環(huán)迭代求解效率低下的問(wèn)題。該方法通過(guò)并行處理能有效避免循環(huán)求解過(guò)程，從而提高計(jì)算效率。

為了實(shí)現(xiàn)上述目的，本發(fā)明的技術(shù)方案是：一種基于并行算法的fetd仿真模擬方法，包括以下步驟：

a.確定需要分析的電真空器件結(jié)構(gòu)；

b.對(duì)步驟a的器件結(jié)構(gòu)進(jìn)行建模，建立對(duì)應(yīng)的幾何結(jié)構(gòu)模型；

c.確定電真空器件結(jié)構(gòu)的電磁學(xué)邊值問(wèn)題的控制微分方程形式；

d.采用四面體單元網(wǎng)格剖分求解區(qū)域；

e.用插值基函數(shù)將控制微分方程中的待求未知量進(jìn)行空間離散展開，并運(yùn)用標(biāo)準(zhǔn)變分原理得到邊值問(wèn)題關(guān)于時(shí)間偏微分的有限元方程組；

f.選擇穩(wěn)定的時(shí)間差分格式(如中心差分格式、newmark-β差分格式)對(duì)步驟e中的有限元方程組進(jìn)行時(shí)間離散，得到邊值問(wèn)題的時(shí)間推進(jìn)方程。

g.采用并行算法計(jì)算步驟f中的時(shí)間推進(jìn)方程的迭代求解過(guò)程。

與現(xiàn)有技術(shù)相比，本發(fā)明的有益效果：利用本發(fā)明提出的一種基于fetd的時(shí)間推進(jìn)方程迭代求解過(guò)程的并行實(shí)現(xiàn)算法，在保證計(jì)算精度的同時(shí)能夠有效提高該過(guò)程計(jì)算效率。

附圖說(shuō)明

圖1是矩形波導(dǎo)結(jié)構(gòu)網(wǎng)格離散后的示意圖。

圖2是本發(fā)明基于并行算法的fetd仿真模擬方法的流程圖。

圖3是fetd的時(shí)間推進(jìn)方程迭代求解過(guò)程串行實(shí)現(xiàn)流程圖。

圖4是fetd的時(shí)間推進(jìn)方程迭代求解過(guò)程并行實(shí)現(xiàn)流程圖。

具體實(shí)施方式

下面結(jié)合附圖和具體實(shí)施方式來(lái)詳細(xì)描述本發(fā)明的技術(shù)方案。

參照流程圖2，一種基于并行算法的fetd仿真模擬方法包括以下步驟：

a.確定需要分析的電真空器件結(jié)構(gòu)。

選取需要分析的電真空器件結(jié)構(gòu)，如波導(dǎo)、矩形窗、盒形窗、同軸窗等，本實(shí)施例選取矩形波導(dǎo)。

b.對(duì)步驟a的器件結(jié)構(gòu)進(jìn)行建模，建立對(duì)應(yīng)的幾何結(jié)構(gòu)模型；

本發(fā)明基于矩形波導(dǎo)結(jié)構(gòu)而言，網(wǎng)格離散后的幾何結(jié)構(gòu)模型如圖1所示。

c.確定電真空器件結(jié)構(gòu)的電磁學(xué)邊值問(wèn)題的控制微分方程形式

對(duì)于矩形波導(dǎo)結(jié)構(gòu)而言，對(duì)該結(jié)構(gòu)電磁傳輸特性的分析是在一定的空間內(nèi)求解maxwell方程組，該問(wèn)題的求解域?yàn)椴▽?dǎo)壁界定的空間。從maxwell方程組出發(fā)，在無(wú)外加電流激勵(lì)的情況下可以得到電場(chǎng)e滿足的微分方程為：

其中，μ為磁導(dǎo)率，ε為介電常數(shù)，σ為電導(dǎo)率，t為時(shí)間。為了使公式(1)中的方程解唯一，需要確定波導(dǎo)入射端口面和出射端口面上的邊界條件。一般邊界條件滿足：

上式中，表示端口面上的單位外法矢量，γ為邊界等效波導(dǎo)納，為已知量。

d.采用四面體單元網(wǎng)格剖分求解區(qū)域。

采用四面體網(wǎng)格剖分求解域是時(shí)域有限元方法中的一種公知過(guò)程，因此本步驟不再詳細(xì)描述。剖分后的求解域被人為分割為多個(gè)三維四面體網(wǎng)格，從而將連續(xù)的幾何結(jié)構(gòu)空間轉(zhuǎn)化為離散的網(wǎng)格空間。

e.用插值基函數(shù)將控制微分方程中的待求未知量進(jìn)行空間離散展開，并運(yùn)用標(biāo)準(zhǔn)變分原理得到邊值問(wèn)題關(guān)于時(shí)間偏微分的有限元方程組。

如公式(1)所示，對(duì)于步驟c中的控制微分方程的求解需要先用插值基函數(shù)將待求未知量進(jìn)行空間離散展開。這個(gè)也是公知的過(guò)程，因此本步驟不再詳細(xì)描述，直接給出通過(guò)標(biāo)準(zhǔn)變分原理將控制微分方程進(jìn)行空間離散后的具體形式：

其中{u}為電場(chǎng)按基函數(shù)展開的展開系數(shù)。各系數(shù)矩陣[t]、[tσ]、[s]、{f}中的矩陣元素滿足：

公式(4)中，v是體積，s是邊界面面積，ni和nj均為插值基函數(shù)，插值基函數(shù)為矢量棱邊基函數(shù)時(shí)i的取值為：i＝1,2,3,4,5,6。同樣j的取值范圍和i一致。

f.選擇穩(wěn)定的時(shí)間差分格式對(duì)步驟e中的有限元方程組進(jìn)行時(shí)間離散，得到邊值問(wèn)題的時(shí)間推進(jìn)方程。

在步驟e的基礎(chǔ)上，對(duì)公式(3)中的方程進(jìn)行時(shí)間離散。這里選取中心差分格式進(jìn)行離散得到時(shí)間推進(jìn)方程為：

其中，δt為時(shí)間步長(zhǎng)，{f}ⁿ為激勵(lì)，{u}ⁿ為第n時(shí)刻的待求未知量。

g.通過(guò)并行算法實(shí)現(xiàn)步驟f中的時(shí)間推進(jìn)方程的迭代求解過(guò)程。

對(duì)于公式(5)中的時(shí)間推進(jìn)方程，在給出激勵(lì)形式、初始值{u}⁰和{u}¹后，可以求解出{u}²的值；再通過(guò){u}¹和{u}²求出{u}³的值，依次類推可以求出后續(xù)任意時(shí)刻的{u}值。這種方式可以看成一個(gè)迭代遞推的過(guò)程，目前處理該過(guò)程采用的方法都是用循環(huán)的方式串行實(shí)現(xiàn)，其流程如圖3所示。流程圖3展示了待求問(wèn)題循環(huán)迭代求解的串行過(guò)程，初始零時(shí)刻時(shí)，{u}⁰,{u}¹,f⁰已知，在滿足循環(huán)條件的情況下進(jìn)入循環(huán)體可以計(jì)算出{u}²的值，接著返回循環(huán)條件繼續(xù)判斷是否滿足該條件，如果滿足循環(huán)條件，則繼續(xù)進(jìn)入循環(huán)體計(jì)算{u}³,{u}⁴，…，{u}ⁿ。循環(huán)條件是通過(guò)迭代次數(shù)進(jìn)行約束，在循環(huán)之前需要知道該約束條件。上述過(guò)程需要通過(guò)不斷的循環(huán)迭代求解從而得到后續(xù)時(shí)刻的值，這種串行計(jì)算的思想在計(jì)算效率上并不能達(dá)到最優(yōu)。

本發(fā)明中步驟a-f為常規(guī)計(jì)算過(guò)程，核心發(fā)明點(diǎn)為采用并行算法替代現(xiàn)有的串行算法，具體實(shí)施過(guò)程為：

由公式(5)，fetd的時(shí)間推進(jìn)方程可以簡(jiǎn)化為：

[a]{u}ⁿ⁺¹＝[b]{u}ⁿ+[c]{u}^n-1+{f}ⁿ(6)

其中[a],[b],[c]為系數(shù)矩陣，并且由公式(5)中的各系數(shù)矩陣得到。這里對(duì)公式(6)兩邊同乘[a]^-1得到：

{u}ⁿ⁺¹＝[a]^-1[b]{u}ⁿ+[a]^-1[c]{u}^n-1+[a]^-1{f}ⁿ(7)

在主模激勵(lì)的情況下，{f}ⁿ滿足：

聯(lián)立公式(7)和公式(8)可得：

將公式(9)進(jìn)行簡(jiǎn)化得到：

{u}ⁿ⁺¹＝[m]n{u}ⁿ(10)

其中：

根據(jù)公式(10)并令[m]¹＝[m]n，[m]²＝[m]n[m]n-1，…，[m]ⁿ＝[m]n[m]n-1…[m]1，則有：

{u}ⁿ＝[m]¹{u}^n-1＝[m]²{u}^n-2＝…＝[m]^n-1{u}¹(11)

也即對(duì)公式(3)的求解轉(zhuǎn)化為公式(11)的解。對(duì)于公式(11)的迭代求解過(guò)程，采用并行算法實(shí)現(xiàn)，其流程圖如圖4所示。流程圖4展示了待求問(wèn)題的并行實(shí)現(xiàn)過(guò)程，在給定初值的情況下，程序?qū)崿F(xiàn)時(shí)需要?jiǎng)?chuàng)建多個(gè)線程，線程的個(gè)數(shù)根據(jù)需要迭代的次數(shù)確定。在計(jì)算過(guò)程中，讓多個(gè)線程同時(shí)工作，例如：線程1計(jì)算{u}²的值，線程2計(jì)算{u}³的值，線程3計(jì)算{u}⁴的值等等。從流程圖4可以看到，各個(gè)線程的工作相互獨(dú)立且同時(shí)進(jìn)行。從加快計(jì)算速度和提高計(jì)算效率的角度考慮，上述這種并行的處理方式能夠解決目前傳統(tǒng)方法中存在的問(wèn)題。這種并行處理迭代求解過(guò)程的方式和串行循環(huán)處理相比在計(jì)算效率上得到優(yōu)化，能夠充分利用計(jì)算機(jī)的資源。

還需要注意的是對(duì)于n×n維矩陣[m]，當(dāng)n取20000時(shí)用數(shù)組的方式存儲(chǔ)該矩陣需要8×20000×20000個(gè)字節(jié)(3gb)的內(nèi)存。這種直接存儲(chǔ)矩陣元素的方式對(duì)內(nèi)存的消耗過(guò)大，在數(shù)值計(jì)算編程過(guò)程中是要避免的。由于矩陣[m]對(duì)稱正定且矩陣中有很多零元素，而這些零元素是不直接參與計(jì)算的，對(duì)計(jì)算產(chǎn)生作用的只有那些非零元素。從考慮節(jié)約內(nèi)存的角度出發(fā)，我們采用稀疏存儲(chǔ)的方式來(lái)存儲(chǔ)矩陣中的非零元素。稀疏存儲(chǔ)的方式主要有按行稀疏存儲(chǔ)和按列稀疏存儲(chǔ)。這里采用按列稀疏存儲(chǔ)的方式來(lái)進(jìn)行，對(duì)于n階m個(gè)非零元素的矩陣來(lái)說(shuō)，用三個(gè)一維數(shù)組*ptr，*ind和*val來(lái)稀疏存儲(chǔ)，ind，val的長(zhǎng)度是m，按行順序記錄每個(gè)非零元的列標(biāo)和數(shù)值；ptr的長(zhǎng)度是n+1，ptr[i]記錄第i行第一個(gè)非零元的位置。最后一個(gè)元素ptr[n]＝m。例如對(duì)于三維矩陣：

則：

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}；

因?yàn)閇m]是按列稀疏方式存儲(chǔ)，在計(jì)算公式(11)時(shí)，可以利用英特爾數(shù)學(xué)核心函數(shù)庫(kù)(intelmkl)里的相關(guān)庫(kù)函數(shù)進(jìn)行求解計(jì)算。用稀疏的方式存儲(chǔ)矩陣并結(jié)合intelmkl里的庫(kù)函數(shù)對(duì)公式(11)中的方程進(jìn)行求解，不僅節(jié)約了大量?jī)?nèi)存空間，更提高了求解效率。