[發(fā)明專利]基于完全矩陣計(jì)算的離散元鄰居搜索及求解方法和系統(tǒng)

在審

• 申請?zhí)枺?/td>: 201911066584.1
• 申請日: 2019-11-04
• 公開（公告）號: CN110781448A
• 公開（公告）日: 2020-02-11
• 發(fā)明（設(shè)計(jì)）人: 劉春;劉輝;施斌;王寶軍;樂天呈;張丹;唐朝生;朱鴻鵠;顧凱
• 申請（專利權(quán)）人: 南京大學(xué)
• 主分類號: G06F17/16
• 分類號: G06F17/16;G06F16/901
• 代理公司: 32249 南京瑞弘專利商標(biāo)事務(wù)所(普通合伙)
• 代理人: 彭雄

本發(fā)明公開了一種基于完全矩陣計(jì)算的離散元鄰居搜索及求解方法和系統(tǒng)，將空間劃分為正交網(wǎng)格元胞；通過完全矩陣運(yùn)算，基于網(wǎng)格元胞的鄰域元胞建立初步接觸搜索矩陣，并通過距離判別篩選接觸信息及虛單元填充方法，得到接觸搜索的鄰居矩陣表示；基于鄰居矩陣及相應(yīng)布爾矩陣實(shí)現(xiàn)離散元程序的求解計(jì)算。本發(fā)明摒棄了循環(huán)語句結(jié)構(gòu)，實(shí)現(xiàn)了離散元鄰居搜索與求解計(jì)算的純矩陣運(yùn)算，進(jìn)而通過GPU并行加速技術(shù)大幅提高三維離散元計(jì)算效率。

技術(shù)領(lǐng)域

本發(fā)明涉及一種鄰居搜索及求解計(jì)算的純矩陣實(shí)現(xiàn)方法，尤其適用于顆粒離散元GPU并行加速的高效計(jì)算方法。

背景技術(shù)

自1979年Cundall提出離散元法以來，離散元法因其天然的非連續(xù)特性被廣泛應(yīng)用到顆粒物質(zhì)材料及大變形問題。離散元法基于牛頓第二定律，通過時(shí)間步迭代算法顯示求解，因此隨著計(jì)算的空間尺度和時(shí)間尺度的增加，計(jì)算單元數(shù)和時(shí)間步迭代次數(shù)逐漸增加，離散元法的計(jì)算量也顯著增加。其中三維接觸搜索與判別是顆粒離散元計(jì)算中最關(guān)鍵的問題之一，其計(jì)算量占比可達(dá)到總計(jì)算量的50％到80％。特別是對于大變形問題(如滑坡)、動力過程，如何在數(shù)百萬單元運(yùn)動過程中，實(shí)時(shí)高效地實(shí)現(xiàn)三維接觸搜索判別一直是制約離散元法在工程尺度推廣與應(yīng)用的瓶頸。

近年來基于GPU的高性能計(jì)算逐漸興起，而傳統(tǒng)的離散元程序往往以串行代碼為基礎(chǔ)，以鏈表的形式存儲接觸信息，不利于代碼的GPU并行計(jì)算。

發(fā)明內(nèi)容

發(fā)明目的：為了克服現(xiàn)有技術(shù)中存在的不足，本發(fā)明提供一種基于完全矩陣計(jì)算的離散元鄰居搜索及求解方法和系統(tǒng)，本發(fā)明通過虛單元填充方法以鄰居矩陣的形式存儲接觸信息，實(shí)現(xiàn)了顆粒離散元程序的純矩陣運(yùn)算，便于程序編寫與GPU并行加速，實(shí)現(xiàn)了高效的離散元鄰居搜索及求解計(jì)算。

技術(shù)方案：為實(shí)現(xiàn)上述目的，本發(fā)明采用的技術(shù)方案為：

一種基于完全矩陣計(jì)算的離散元鄰居搜索方法，包括以下步驟：

步驟1，單元編號及信息存儲：將單元順序編號，記從1到m，同時(shí)添加虛單元，其編號m+1，通過編號唯一標(biāo)識單元；通過虛單元填充得到接觸搜索的鄰居矩陣表示，鄰居矩陣為m行，每一行記錄單元可能接觸的單元編號，長短不一部分用虛單元編號填充，單元信息存儲為單元屬性數(shù)組，記為aX、aY、aZ和aR，其中，編號m+1單元為虛單元，aX、aY、aZ為三維坐標(biāo)，aR為半徑；

步驟2，網(wǎng)格劃分及編號：

步驟21，輸入不包括虛單元的單元三維坐標(biāo)及半徑數(shù)組，記三維坐標(biāo)數(shù)組為X,Y,Z，半徑數(shù)組為R，其均為m行的一維數(shù)組；

步驟22，根據(jù)三維坐標(biāo)數(shù)組確定網(wǎng)格邊界及邊長：

p方向網(wǎng)格兩側(cè)邊界：

Pmin＝min(P)-minR/2

Pmax＝max(P)+minR/2

其中，[Pmin，Pmax]為p方向兩側(cè)邊界，[min(P),max(P)]為p方向坐標(biāo)數(shù)組最小最大值，minR為最小單元半徑，P表示x,y,z方向；

記[Xmin,Xmax],[Ymin,Ymax],[Zmin,Zmax]為x,y,z方向網(wǎng)格兩側(cè)邊界；

網(wǎng)格元胞邊長：

gSide＝2*maxR+dSide

其中，gSide為網(wǎng)格元胞邊長，maxR為單元最大半徑，dSide為可調(diào)值；

以最小值邊界為起點(diǎn)，gSide為網(wǎng)格元胞邊長，劃分空間網(wǎng)格；則p方向元胞數(shù)PNum：

其中，ceil表示向上取整運(yùn)算，元胞沿p方向順序編號，記從1到PNum；