[發(fā)明專利]體和面混合GPU并行的計算電磁學(xué)DGTD方法

有效

• 申請?zhí)枺?/td>: 201811476653.1
• 申請日: 2018-12-05
• 公開（公告）號: CN109726441B
• 公開（公告）日: 2022-05-03
• 發(fā)明（設(shè)計）人: 徐立;王晗;楊中海;李斌
• 申請（專利權(quán)）人: 電子科技大學(xué)
• 主分類號: G06F30/23
• 分類號: G06F30/23;G06T17/20;G06T1/20;G06F111/10
• 代理公司: 電子科技大學(xué)專利中心 51203
• 代理人: 閆樹平

1.體和面混合GPU并行的計算電磁學(xué)DGTD方法，具體步驟如下：

S1、在主機(jī)端為讀取的三角形面網(wǎng)格數(shù)組按照邊界條件類型進(jìn)行排序，同一邊界條件的面的編號相連，內(nèi)部面置于最后，之后計算質(zhì)量矩陣的逆的系數(shù)矩陣，大小為d_i*d_i，計算面通量矩陣的系數(shù)矩陣，大小為d_f*d_f，d_i為每單元的DOF，d_f為每個面的DOF；并對電場數(shù)組:E_x、E_y、E_z和磁場數(shù)組H_x、H_y、H_z進(jìn)行初始化，每個場的數(shù)組大小為d_i*K，K為單元總數(shù)；其中DOF為自由度；

S2、為GPU端的電磁場數(shù)組，所有單元的體右端項，所有面的面右端項以及電磁場數(shù)組分配好GPU全局內(nèi)存空間，為質(zhì)量矩陣的逆的系數(shù)矩陣、面通量矩陣的系數(shù)矩陣分配好GPU常量內(nèi)存，內(nèi)存分配具體數(shù)量為：

體右端項三個分量的數(shù)組rhsx、rhsy、rhsz的大小分別為d_i*K，面右端項又分為一個面在其左側(cè)單元的右端項的三個分量rhsx_l、rhsy_l、rhsz_l和在其右側(cè)單元的右端項的三個分量rhsx_r、rhsy_r、rhsz_r，每個數(shù)組的大小為d_i*K_f；質(zhì)量矩陣的逆矩陣和通量矩陣的系數(shù)矩陣大小分別為d_i*d_i，d_f*d_f；K_f為面的總數(shù)，d_f為每個面的DOF，并將在主機(jī)端初始化完畢的電場、磁場數(shù)組從主機(jī)內(nèi)存復(fù)制到GPU顯存中；

S3、在主機(jī)端判斷時間步n是否達(dá)到仿真步數(shù)上限，若是，結(jié)束計算，釋放所有已分配的GPU顯存空間，否則繼續(xù)；

S4、在GPU端，先計算每個單元的剛度矩陣，每個線程需要計算剛度矩陣的一行，大小為d_i，再計算每個單元的電場體右端項為第i單元的剛度矩陣，為第i單元按照基函數(shù)展開的磁場矢量的展開系數(shù)，且線程的分配方式如下：

依據(jù)每個線程計算一個DOF的原則分配，每個線程塊負(fù)責(zé)處理N_T個單元，使得一個線程塊內(nèi)的線程數(shù)量取：

上式中16為CUDA線程束大小的一半，符號代表A向下取整，本步驟計算結(jié)果存儲在對應(yīng)的體右端項三個坐標(biāo)分量的數(shù)組rhsx、rhsy、rhsz中；

S5、在GPU端，先由通量矩陣的系數(shù)矩陣計算出每個面的通量矩陣，再計算每個面對其左側(cè)單元、右側(cè)單元的電場面右端項：和v∈{x,y,z}，為第i單元的通量矩陣，其中L表示面的左側(cè)單元，R表示面的右側(cè)單元，n為單位外法向量，為數(shù)值通量；且線程的分配方式如下：

依據(jù)每個線程計算一個DOF的原則分配，每個線程塊負(fù)責(zé)處理N_F個面的左、右右端項，使得一個線程塊內(nèi)的線程數(shù)量?。?/p>

本步驟的計算結(jié)果根據(jù)是對左側(cè)還是右側(cè)單元的貢獻(xiàn)而分別存儲在左、右面對應(yīng)的右端項三個坐標(biāo)分量的數(shù)組rhsx_l、rhsy_l、rhsz_l和rhsx_r、rhsy_r、rhsz_r中；

S6、在GPU端，先由質(zhì)量矩陣的逆的系數(shù)矩陣計算每個單元質(zhì)量矩陣的逆，再組合每個單元電場的體右端項和面右端項為總右端項，并用質(zhì)量矩陣的逆乘以總右端項得到為第i單元電場的質(zhì)量矩陣，為第i單元的第f個面的通量矩陣，n_if為第i單元的第f個面的單位外法向量，為此面的數(shù)值通量；最后計算得到新時刻的電場；

本步驟核函數(shù)的線程分配與S4一致，計算完畢后將GPU端的電場數(shù)組使用cudaMemcpyAsync()函數(shù)異步拷貝到主機(jī)端進(jìn)行輸出或者后處理；

S7、在GPU端，先計算每個單元的剛度矩陣，每個線程需要計算剛度矩陣的一行，再計算每個單元的磁場體右端項為第i單元按照基函數(shù)展開的電場矢量的展開系數(shù)，且線程的分配方式同S4，本步驟核函數(shù)計算結(jié)果存儲在對應(yīng)的體右端項三個坐標(biāo)分量的數(shù)組rhsx、rhsy、rhsz中；

S8、在GPU端，先由通量矩陣的系數(shù)矩陣計算每個面的通量矩陣，再計算每個面對其左側(cè)單元、右側(cè)單元的磁場面右端項：和v∈{x,y,z}，且線程的分配方式同S5，邊界條件面的處理同S5，本步驟核函數(shù)的計算結(jié)果根據(jù)是對左側(cè)還是右側(cè)單元的貢獻(xiàn)而分別存儲在對應(yīng)的左、右面的右端項三個坐標(biāo)分量的數(shù)組rhsx_l、rhsy_l、rhsz_l和rhsx_r、rhsy_r、rhsz_r中；

S9、在GPU端，先由質(zhì)量矩陣的逆矩陣的系數(shù)矩陣計算每個單元質(zhì)量矩陣的逆，再組合每個單元磁場的體右端項和面右端項為總右端項，并用剛度矩陣的逆乘總右端項得到為第i單元磁場的質(zhì)量矩陣，為此面的數(shù)值通量，最后計算得到新時刻的磁場；本步驟核函數(shù)的線程分配方式同S4，計算完畢后將GPU端的磁場數(shù)組使用cudaMemcpyAsync()異步拷貝到主機(jī)端進(jìn)行輸出或者后處理，一個時間步長計算完成，轉(zhuǎn)S3；

所述步驟S4和S7在計算體右端項時，每個線程從全局內(nèi)存讀取一個DOF的場；

所述步驟S5和S8的線程分配是基于面網(wǎng)格的，每個面對其左右兩個單元的右端項都有不同的貢獻(xiàn)，如果用下標(biāo)L和R來分別標(biāo)識一個面的左、右兩個單元的物理量，n_LR表示當(dāng)前計算的面上由左單元指向右單元的單位法向量，而n_RL＝-n_LR表示當(dāng)前面上由右單元指向左單元的單位法向量，則一個面關(guān)于的迎風(fēng)通量可以寫為

進(jìn)而面右端項可以寫為

將n_RL＝-n_LR帶入公式(4)可以得到

其中Z為阻抗，S_LL,S_RR,S_LR,S_RL分別代表由大小為d_f*d_f擴(kuò)展為d_i*d_i的面通量矩陣，擴(kuò)展方式為根據(jù)該面的在左右單元中的局部編號，將原通量矩陣的每行每列放到擴(kuò)展矩陣中該面的自由度在兩個單元中所對應(yīng)的自由度的行列中，下標(biāo)代表矩陣是相對于左右哪個單元擴(kuò)展的，第一個下標(biāo)是對行的擴(kuò)展，第二個下標(biāo)是對列的擴(kuò)展；

在這其中，每個線程從全局內(nèi)存中讀取一個面通量矩陣的一行，每個線程先讀取左單元中一個DOF的電磁場，在計算完公式(5)中需要使用左單元場的四項之后，讀取右單元中一個DOF的電磁場并計算其余四項。