[發明專利]一種基于GPU并行的探地雷達電磁波數值模擬計算方法

有效

• 申請號: 202010182487.5
• 申請日: 2020-03-16
• 公開（公告）號: CN111324972B
• 公開（公告）日: 2023-02-24
• 發明（設計）人: 方宏遠;雷建偉;薛冰寒;李銀萍;胡浩幫;張超
• 申請（專利權）人: 鄭州大學
• 主分類號: G06F30/20
• 分類號: G06F30/20;G06F111/10
• 代理公司: 鄭州銀河專利代理有限公司 41158
• 代理人: 陳玄

1.一種基于GPU并行的探地雷達電磁波數值模擬計算方法，其特征在于，包括如下步驟：

步驟1、將待求解的地下結構模型區域劃分為大小相同的網格，電場和磁場定義在網格正中心；在計算機的CPU上分配存儲空間，用于儲存二維平面地下結構離散區域的長度和寬度、離散區域的網格步長和時間步長、磁場的磁導率、電場的電導率、介電常數參數，并對所述參數初始化；同時在GPU上也分配存儲空間用于存儲：二維平面地下結構離散區域的長度和寬度、離散區域的網格步長和時間步長、磁場的磁導率、電場的電導率、介電常數參數，并對這些參數初始化；其中在CPU上定義地下結構模型區域的離散劃分的各網格磁場和電場的結構體數組，用于在GPU上存儲通過磁場和電場公式(1)并行計算得到的各時間步各網格的磁場值和電場值，所述的時間步，就是把整體的時間分為n段，每一段就是一個時間步，在第一個時間步內，計算各網格的磁場和電場，然后保存數據，在第二個時間步內再全部計算各網格的磁場和電場，在計算第二個時間步時需要第一個時間步的磁場和電場的值，依次一直計算到第n個時間步；這些并行計算得到的場值通過傳輸函數cudaMemcpy傳輸到CPU上定義的各磁場和電場的結構體數組，其中所述磁場和電場的初始值都為0；

步驟2、利用下面公式計算地下結構模型區域離散各網格的磁場和電場的值，二維辛計算方法的磁場A和電場U的計算公式如式(1)所示：

式中，σ表示的電導率，ε表示介電常數，μ表示磁導率，Δt為模擬計算的時間步長，Δx為二維地下結構模型區域離散網格在x方向的長度，Δy為二維地下結構模型區域離散網格在y方向的長度，i為二維平面地下結構模型模擬區域離散網格在y方向的編號，j為二維平面地下結構模型模擬區域離散網格在x方向的編號，表示網格(i，j)處磁場A第n個時間步長時的值，表示網格(i，j)處電場U第n個時間步長時的值，表示網格(i，j)處磁場A第n+1個時間步長時的值，表示網格(i，j)處電場U第n+1個時間步長時的值；

步驟3、在模擬區域網格上加載探地雷達脈沖，然后使用GPU上設置的線程并行計算公式(1)中第一式的磁場A，得到各網格處的磁場A的值；

步驟4、利用步驟3計算得到的磁場A的值，在GPU上利用設置的并行線程計算公式(1)第二式中的電場U，得到各網格處的電場U的值；最后進行電場U的邊界處理，采用坐標伸縮完全匹配層吸收邊界條件，二維辛計算方法的磁場U的坐標伸縮完全匹配層吸收邊界條件的計算公式如式(2)所示：

式中κ_x為x方向的有效延伸因子，κ_y為y方向的有效延伸因子，其中CA和CB為公式系數，具體為

和為中間變量函數可由如下公式計算

式中的a_x為x方向第一公式系數，b_x為x方向第一公式系數，a_y為y方向第一公式系數，b_y為y方向第二公式系數，系數具體求解式如下：

式中，式中κ_i為場值在i方向的有效延伸因子，κ_max為κ_i的最大值，κ_max取5；α_i為場值在i方向的坐標延伸因子的自由度，α_max為α_i的最大值，α_max取0.008；σ_i為坐標伸縮完全匹配層吸收邊界區域中在i方向的電導率，σ_max為σ_i的最大值，其計算公式為式中的m＝4，δ為離散網格的長度；d為坐標伸縮完全匹配層吸收層厚度，d^m是d的m次方；

步驟5、將計算網格的電場和磁場的值存儲到GPU全局內存，并重復步驟3和步驟4，計算二維地下結構模型區域離散每個網格不同時間步的電場A和磁場U，從時間步T＝0迭代直到最大時間步數T＝Tmax，Tmax為模擬計算的最大時間步數；

步驟6、將步驟5計算得到的各網格的電場和磁場的值從GPU設備上拷貝到CPU上，對電場和磁場的值進行存儲，將電場的值導入到MATLAB中，通過MATLAB畫圖分析探地雷達電磁波在地下結構中的傳播規律。