許衛(wèi)明， 王建華

(1.馬鞍山師范高等?？茖W(xué)校，安徽 馬鞍山 243041; 2.菏澤市城市管理局,山東 菏澤 274000)

?

OpenCL離散元顆粒系統(tǒng)的優(yōu)化設(shè)計(jì)*

許衛(wèi)明1， 王建華2

(1.馬鞍山師范高等?？茖W(xué)校，安徽 馬鞍山 243041; 2.菏澤市城市管理局,山東 菏澤 274000)

針對傳統(tǒng)離散元顆粒系統(tǒng)進(jìn)行大量顆粒運(yùn)算時(shí)，出現(xiàn)的運(yùn)算時(shí)間長、運(yùn)算效果不理想的問題，通過OpenCL算法對其進(jìn)行優(yōu)化改進(jìn).將大量顆粒運(yùn)算分為顆粒間碰撞運(yùn)算和顆粒與邊界的碰撞運(yùn)算.顆粒間碰撞的運(yùn)算利用OpenCL編程實(shí)現(xiàn)，與邊界碰撞的運(yùn)算利用CPU計(jì)算實(shí)現(xiàn)，同時(shí)在OpenCL算法運(yùn)行時(shí)，根據(jù)不同的數(shù)據(jù)類型，劃分不同的存儲器，在內(nèi)存訪問上提高算法的運(yùn)行速度.通過測試，優(yōu)化后的算法運(yùn)算效率明顯提升.

顆粒系統(tǒng)；OpenCL；GPU；內(nèi)存訪問

離散元方法的基本思想是把不連續(xù)體分離成剛性元素集合，利用運(yùn)動方程來描述剛性元素運(yùn)動情況，通過迭代的方法對運(yùn)動方程進(jìn)行求解，從而得出不連續(xù)體的運(yùn)動形態(tài).最初該方法是針對非連續(xù)介質(zhì)的力學(xué)行為的[1].在離散元的顆粒系統(tǒng)中，由于顆粒數(shù)量多，存在著大量的顆粒間的碰撞和顆粒與邊界的碰撞，對系統(tǒng)的力學(xué)進(jìn)行計(jì)算，需處理、跟蹤大量的信息，僅憑CPU的運(yùn)算，效果太差.近年來隨著GPU科學(xué)計(jì)算能力的逐步提高，尤其是GPU的并行計(jì)算能力的發(fā)展，將計(jì)算機(jī)的計(jì)算能力大幅提升，使得在對顆粒系統(tǒng)中數(shù)億的顆粒模擬計(jì)算成為可能，但整體運(yùn)算時(shí)間還是較長，本文對基于OpenCL離散元顆粒系統(tǒng)進(jìn)行研究，通過在算法上適應(yīng)GPU構(gòu)架，同時(shí)在內(nèi)存訪問上進(jìn)行優(yōu)化，提升離散元的仿真效率.

1　OpenCL簡介

OpenCL是由蘋果公司和AMD、NVIDIA、Intel等公司制定的一個(gè)工業(yè)標(biāo)準(zhǔn)和規(guī)范，是面向CPU、GPU與其他處理器組合構(gòu)成的計(jì)算機(jī)編程的行業(yè)標(biāo)準(zhǔn).該標(biāo)準(zhǔn)可通過公布硬件來提高編程的可移植性.平臺模型如圖1所示，平臺包括各種OpenCL設(shè)備及一個(gè)宿主機(jī)，宿主機(jī)與各個(gè)OpenCL設(shè)備之間的交互通信，通過宿主機(jī)定義內(nèi)核，建立邏輯關(guān)系，對各個(gè)設(shè)備下達(dá)命令的方式來完成[2].

圖1　OpenCL平臺模型

OpenCL定義了五種內(nèi)存區(qū)域，首先為宿主機(jī)內(nèi)存，該內(nèi)存只有宿主機(jī)可見；全局內(nèi)存，所有設(shè)備可見，并允許讀寫操作；常量內(nèi)存，用于存儲常量，只允許讀；局部內(nèi)存，存儲不同設(shè)備的公用常量；私有內(nèi)存，計(jì)算中一個(gè)工作組的專用內(nèi)存.在利用OpenCL進(jìn)行編程時(shí)一般遵循以下步驟：1)查支持OpenCL的硬件設(shè)備的參數(shù)、特性；2)創(chuàng)建OpenCL平臺的內(nèi)核函數(shù)；3)在設(shè)備上創(chuàng)建上下邏輯關(guān)系，管理程序執(zhí)行時(shí)的內(nèi)核對象；4)傳遞所需的參數(shù)；5)執(zhí)行內(nèi)核函數(shù).

2　基于 OpenCL 離散元顆粒系統(tǒng)仿真優(yōu)化設(shè)計(jì)

離散元顆粒系統(tǒng)仿真是將顆粒系統(tǒng)中的顆粒假象為獨(dú)立的球形或橢球型顆粒，計(jì)算這些顆粒之間碰撞的受力情況，計(jì)算位移與速度，得出下一時(shí)步的位置和速度，循環(huán)此過程直到計(jì)算所有顆粒停止運(yùn)動為止，或到達(dá)規(guī)定運(yùn)算時(shí)步，仿真結(jié)束.本次的優(yōu)化運(yùn)算首先將顆粒的碰撞分為兩類：顆粒與顆粒之間的碰撞、顆粒與邊界之間的碰撞.

2.1顆粒之間的碰撞計(jì)算

顆粒間的碰撞計(jì)算流程如圖2所示.在OpenCL平臺上創(chuàng)建上下邏輯關(guān)系，選擇硬件設(shè)備，創(chuàng)建命令隊(duì)列，并對命令在隊(duì)列中排隊(duì)，從.cl文件中加載OpenCL內(nèi)核代碼.在內(nèi)存中為內(nèi)核函數(shù)的參數(shù)分配空間[3].對于顆粒系統(tǒng)，函數(shù)的主要參數(shù)為：力學(xué)常量，顆粒的數(shù)量、坐標(biāo)、大小、速度和角速度等信息.由于內(nèi)核函數(shù)不能命名指針，只能將顆粒的信息轉(zhuǎn)換成一維數(shù)組，核函數(shù)最終得出的結(jié)果為該時(shí)步下的顆粒碰撞受力情況，然后通過分析顆粒的受力情況來計(jì)算顆粒的運(yùn)動狀態(tài)，得出顆粒的速度，位置信息等狀態(tài)的變化.

圖2　顆粒間碰撞計(jì)算流程

2.1.1接觸判斷與位移增量計(jì)算

首先建立局部坐標(biāo)系，以遍歷球顆粒i的球心為坐標(biāo)原點(diǎn)，與顆粒j的球心連線作為x軸，選取x軸的垂直平面與X-Y平面平行的直線為y軸，再由右手螺旋法則確定z軸，建立空間坐標(biāo)系，如圖3所示.

圖3　顆粒間局部坐標(biāo)系

全局坐標(biāo)與局部坐標(biāo)的轉(zhuǎn)換矩陣為

在轉(zhuǎn)換矩陣中Lij表示兩球i和j的球心距，(Xi,Yi,Zi)、(Xj,Yj,Zj)表示兩球在全局坐標(biāo)系中的坐標(biāo)，(l，m，n)是x軸在全局坐標(biāo)內(nèi)的方向向量.進(jìn)行接觸判斷時(shí)，要將該顆粒與比該顆粒編號大的顆粒進(jìn)行判斷，計(jì)算過程通過GPU并行進(jìn)行計(jì)算，判定球體是否接觸的方法為：計(jì)算兩顆球的球心距與兩個(gè)球的半徑和的大小，當(dāng)球心距大于半徑和時(shí)，兩球沒有接觸，當(dāng)球心距小于半徑和時(shí)，兩球接觸.

2.1.2力學(xué)計(jì)算

兩球發(fā)生碰撞時(shí)，根據(jù)受力分析，顆粒的受力包括法向接觸力、切向接觸力和回轉(zhuǎn)力矩.力學(xué)的計(jì)算過程在計(jì)算接觸判斷的線程中完成.

在t時(shí)步下，發(fā)生碰撞的的顆粒i在x軸方向的法向接觸力為

(1)

在公式中kx為剛性系數(shù)，Cij為兩球疊加量，hx為法向阻尼系數(shù).

在t時(shí)步下沿局部坐標(biāo)系的y、z的作用力為切向彈性力與切向阻尼力的和，則沿y軸的切向接觸力為

(2)

根據(jù)相同的原理可得沿z軸的切向接觸力

(3)

在t時(shí)步下，回轉(zhuǎn)力矩為彈性分量和阻尼分量的和，則回轉(zhuǎn)力矩的計(jì)算公式為

M(t)=Mt-Δt-krΔφ-hrΔφ/Δt

(4)

kr,hr分別為回轉(zhuǎn)度系數(shù)和阻尼系數(shù).

以上的接觸判斷及力學(xué)計(jì)算，使用OpenCL算法，在GPU中完成的，選取平臺并創(chuàng)建命令隊(duì)列，主程序通過命令隊(duì)列對核函數(shù)進(jìn)行任務(wù)分配和操作運(yùn)算.GPU的多線程并行運(yùn)算與CPU的不同，每個(gè)線程訪問內(nèi)存都是連續(xù)的.

2.2顆粒邊界的碰撞計(jì)算

顆粒系統(tǒng)的邊界分為好多種，如遠(yuǎn)行平面、扇形平面等，在進(jìn)行顆粒與邊界碰撞時(shí)要區(qū)分對待.

2.2.1接觸判斷

判斷過程中使用虛函數(shù)P2BContact()返回接觸標(biāo)志，以此來判定與邊界是否接觸.判定的方法為比較球心到邊界的距離和球顆粒的半徑，當(dāng)半徑大于到邊界的距離時(shí)，則判定為接觸，否則不接觸.

設(shè)平面方程ax+by+cz+e=0，邊界的半徑為R，則球心到邊界的距離為

(5)

2.2.2顆粒與邊界的力學(xué)計(jì)算

與邊界的力學(xué)計(jì)算同顆粒間的力學(xué)計(jì)算大致相同，主要的區(qū)別在于顆粒相對于邊界自身的質(zhì)量可以忽略不計(jì)，因此邊界不會因?yàn)轭w粒的碰撞而產(chǎn)生運(yùn)動[1].計(jì)算之前要根據(jù)不同的材質(zhì)建立數(shù)據(jù)庫，方便在程序錄入的情況下調(diào)用相對應(yīng)的參數(shù).

(6)

(7)

計(jì)算結(jié)束后對每個(gè)顆粒的作用力和力矩進(jìn)行相加之后再與核函數(shù)返回計(jì)算結(jié)果相加，根據(jù)受力情況計(jì)算顆粒的運(yùn)動狀態(tài).

3　內(nèi)存訪問優(yōu)化

在OpenCL算法運(yùn)行時(shí)，好的內(nèi)存訪問模式可以加快讀、寫緩沖區(qū)的速度，從而提升計(jì)算運(yùn)行速度和效率.傳統(tǒng)的CPU優(yōu)化內(nèi)存的方式不適應(yīng)OpenCL設(shè)備內(nèi)存，必須探尋新的內(nèi)存優(yōu)化技術(shù)，來提高OpenCL設(shè)備的運(yùn)行速度.

3.1全局內(nèi)存訪問

在OpenCL的核函數(shù)中應(yīng)用比較廣泛的為全局內(nèi)存，決定全局內(nèi)存性能的因素有訪存次數(shù)和訪存模式.訪存次數(shù)就是核函數(shù)中全局內(nèi)存被讀、寫的次數(shù)，次數(shù)過度就會降低程序整體的運(yùn)行效率，而在這些內(nèi)存中，被訪問的次數(shù)最多的數(shù)據(jù)一般為工作組的共享數(shù)據(jù)，解決的方法為利用局部內(nèi)存訪問速度比全局內(nèi)存快的特點(diǎn)，將共享數(shù)據(jù)寫入該工作組的局部內(nèi)存中去[5].訪存模式，GPU中一般采用多線程分段處理，CPU的分段遍歷方式為每個(gè)線程都順序訪問數(shù)組的每一個(gè)元素，如果采用這種方式，在訪問第一個(gè)元素時(shí)，破壞了GPU數(shù)據(jù)訪問的局部要求，會對整個(gè)運(yùn)算效率造成影響，會降低運(yùn)算性能，解決的方法是采用合并讀取的方式，這樣可以將數(shù)組按線程分為不同的運(yùn)算區(qū)域，讀取的都是數(shù)組相鄰的元素，這樣既符合了OpenCL的內(nèi)存局部性，又有效的提高了緩沖區(qū)數(shù)據(jù)的讀、寫速度.

3.2局部內(nèi)存訪問

局部內(nèi)存是GPU每個(gè)線程組上都能快速訪問的內(nèi)存區(qū)域，訪問速度要比全局內(nèi)存速度快，利用局部內(nèi)存可以減少對全局內(nèi)存的讀、寫操作.程序運(yùn)行中如果無法將數(shù)據(jù)分段合并時(shí)，將這些數(shù)據(jù)放入局部內(nèi)存中，這樣可以提高緩沖區(qū)的讀寫速度，從而提高運(yùn)行的效率.存儲體沖突時(shí)降低局部內(nèi)存訪存性能下降的主要因素，本文通過采用均勻訪存的方法盡量避免沖突的發(fā)生.

均勻訪存的方法基本思想為，任意兩個(gè)相鄰的線程同時(shí)訪問局部內(nèi)存中位置間隔相同的區(qū)域，核心就是訪問間隔，間隔不同，訪問的區(qū)域就不同，這樣的訪問也會造成不同的線程訪問同一個(gè)內(nèi)存區(qū)域[6].OpenCL對局部內(nèi)存的讀寫映射到GPU的LDS上，OpenCL可以將32組線程為一組并行訪問.間隔為2的線程訪問如圖4所示.

圖4　間隔為2的均勻訪問

從圖4可看出存儲0和2處發(fā)生沖突，如果將間隔設(shè)為3，訪問情況如圖5所示.

圖5　間隔為3的均勻訪問

從圖中可以得出每一個(gè)線程都有不同的存儲體，整個(gè)周期內(nèi)訪問均沒出現(xiàn)存儲體沖突.因此只要選擇的間隔合適，采用均勻訪存的方法可以避免存儲體沖突的發(fā)生，提高整個(gè)程序運(yùn)行的效率.

4　結(jié)論

通過OpenCL算法對傳統(tǒng)離散元顆粒系統(tǒng)進(jìn)行優(yōu)化改進(jìn)，將大量顆粒運(yùn)算分為顆粒間碰撞運(yùn)算和顆粒與邊界的碰撞運(yùn)算，并利用OpenCL編程實(shí)現(xiàn)顆粒間碰撞的運(yùn)算，利用CPU計(jì)算實(shí)現(xiàn)邊界碰撞的運(yùn)算，同時(shí)在OpenCL算法運(yùn)行時(shí)，根據(jù)不同的數(shù)據(jù)類型，劃分不同的存儲器，在內(nèi)存訪問上提高算法的運(yùn)行速度.極大的發(fā)揮了GPU的并行計(jì)算能力，提升離散元顆粒系統(tǒng)的仿真效率.

[1]趙啦啦,趙躍民,劉初升,等.濕顆粒堆力學(xué)特性的離散元法模擬研究[J]. 物理學(xué)報(bào)， 2014(03):1-9.

[2]季順迎,趙金鳳,狄少丞,等.面向環(huán)境力學(xué)的離散元分析軟件研發(fā)和工程應(yīng)用[J]. 計(jì)算機(jī)輔助工程， 2014(1):69-75.

[3] 詹云,趙新燦,譚同德.基于OpenCL的異構(gòu)系統(tǒng)并行編程[J]. 計(jì)算機(jī)工程與設(shè)計(jì), 2012(11):4191-4195.

[4] 陳鋼,吳百鋒.面向OpenCL模型的GPU性能優(yōu)化[J]. 計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào), 2011(04):571-581.

[5] 韓博,周秉鋒.GPGPU性能模型及應(yīng)用實(shí)例分析[J]. 計(jì)算機(jī)輔助設(shè)計(jì)與圖形學(xué)學(xué)報(bào), 2009(09).

[6] 徐泳,孫其誠,張凌,等.顆粒離散元法研究進(jìn)展[J]. 力學(xué)進(jìn)展, 2003(02):251-259.

On Optimization Design of OPENCL Discrete Element Particle System

XU Wei-ming1, WANG Jian-hua2

(1.Maanshan Teacher’s College, Maanshan Anhui 243041, China 2. Heze City Administration Bureau, Heze Shandong 274000, China)

Due to the long operation time and unsatisfactory effect of traditional discrete element particle system carrying out a large number of granular computing, OPENCL algorithm is used for the optimization. The large number of particles are divided into inter particle collision and the collision between the particles and the boundary. The operation of inter particle collision is realized by OPENCL programming and the operation of the boundary collision is realized by CPU. At the same time, different memory is partitioned according to different data types in OPENCL algorithm. By testing, the optimized algorithm is obviously improved.

particle system; OPENCL; GPU; memory access

1673-2103(2016)02-0023-05

2015-12-07

安徽省一般教學(xué)研究項(xiàng)目(2013jyxm300 )

許衛(wèi)明,男(1981-)，安徽懷寧人，講師，碩士，研究方向：計(jì)算機(jī)網(wǎng)絡(luò)

F416.67

A