范小康，夏澤宇，龍思凡，楊燦群

(國(guó)防科技大學(xué)計(jì)算機(jī)學(xué)院,湖南 長(zhǎng)沙 410073)

1 引言

數(shù)值模擬技術(shù)是工程應(yīng)用和科學(xué)研究的一個(gè)重要手段。與實(shí)驗(yàn)手段相比，數(shù)值模擬成本低，并且可以模擬很多難以開(kāi)展實(shí)驗(yàn)或無(wú)法開(kāi)展實(shí)驗(yàn)的情況，因而受到越來(lái)越多研究人員的青睞。光滑粒子流體動(dòng)力學(xué)SPH(Smoothed Particle Hydrodynamics)[1 - 3]方法是一種無(wú)網(wǎng)格的粒子方法，近年來(lái)逐漸成為數(shù)值模擬技術(shù)領(lǐng)域的一個(gè)研究熱點(diǎn)。

傳統(tǒng)的數(shù)值模擬技術(shù)以網(wǎng)格作為基礎(chǔ)，將實(shí)際問(wèn)題轉(zhuǎn)化為數(shù)學(xué)描述的離散形式，即物理控制方程組，然后采用計(jì)算機(jī)求解。比較具有代表性的方法包括有限差分法FDM(Finite Differential Method)[4]、有限體積法FVM(Finite Volume Method)[5]和有限元法FEM(Finite Element Method)[6]，已經(jīng)在計(jì)算流體力學(xué)CFD(Computational Fluid Dynamics)和計(jì)算固體力學(xué)CSM(Computational Solid Mechanics)中取得了廣泛應(yīng)用。傳統(tǒng)基于網(wǎng)格的方法存在著固有的缺陷，復(fù)雜幾何形狀的網(wǎng)格構(gòu)造是一個(gè)難點(diǎn)，此外，基于網(wǎng)格的方法也難以處理大變形問(wèn)題。

近年來(lái)，無(wú)網(wǎng)格方法受到了越來(lái)越多的重視。無(wú)網(wǎng)格方法的主要思想是通過(guò)使用一系列任意分布的節(jié)點(diǎn)或者粒子來(lái)求解具有各種各樣邊界條件的積分方程或偏微分方程組。無(wú)網(wǎng)格方法的優(yōu)點(diǎn)在于免去了網(wǎng)格生成的麻煩，并且避免了網(wǎng)格扭曲與網(wǎng)格重構(gòu)的問(wèn)題。SPH是一種無(wú)網(wǎng)格的拉格朗日粒子方法，該方法的基本思想是將連續(xù)的流體(或固體)用相互作用的質(zhì)點(diǎn)組來(lái)描述，各個(gè)質(zhì)點(diǎn)上承載各種物理量，包括質(zhì)量、位置坐標(biāo)、速度等，通過(guò)求解質(zhì)點(diǎn)組的動(dòng)力學(xué)方程和跟蹤每個(gè)質(zhì)點(diǎn)的運(yùn)動(dòng)軌道，求得整個(gè)系統(tǒng)的力學(xué)行為。SPH方法的粒子特性使得其在處理大變形、運(yùn)動(dòng)物質(zhì)表面以及自由表面等問(wèn)題時(shí)優(yōu)勢(shì)明顯，已被諸多商用軟件采用。

作為一種粒子方法，SPH方法在每一個(gè)迭代步都需要計(jì)算大量的粒子間相互作用，如何提高SPH方法的計(jì)算效率是目前的研究難點(diǎn)，也具有非常大的實(shí)際意義。針對(duì)SPH方法的計(jì)算效率問(wèn)題，國(guó)內(nèi)外也開(kāi)展了大量的研究工作[7 - 11]。作為一種顯式求解的粒子方法，SPH天然地適合采用并行結(jié)構(gòu)求解，目前已涌現(xiàn)出大量基于分布式存儲(chǔ)的高性能計(jì)算平臺(tái)、共享存儲(chǔ)結(jié)構(gòu)的并行SPH計(jì)算方法，有效地提高了SPH方法的計(jì)算效率，使得SPH可處理的問(wèn)題規(guī)模越來(lái)越大。然而，以上方法研究的都是SPH方法中的進(jìn)程級(jí)或線程級(jí)并行，卻鮮有工作探索SPH方法中的數(shù)據(jù)級(jí)并行。

近年來(lái)ARM面向高性能計(jì)算推出了可伸縮矢量擴(kuò)展SVE(Scalable Vector Extension)[12]。SVE并不是ARM NEON的簡(jiǎn)單擴(kuò)展，而是基于ARM AArch64架構(gòu)的下一代SIMD指令集。與NEON相比，SVE可以獲得超過(guò)7倍的加速比，極大地提升了性能[12]。SVE允許硬件設(shè)計(jì)師根據(jù)應(yīng)用需求定義向量寄存器寬度，最大可支持2 048位的向量寬度。SVE支持VLA(Vector Length Agnostic)編程方式，編程人員或編譯器無(wú)需在程序中指定向量寬度，程序被編譯成SVE指令后，可在不同向量寬度的CPU上兼容運(yùn)行，而無(wú)需針對(duì)不同向量寬度的CPU重新編譯。SVE在高性能計(jì)算領(lǐng)域擁有非常好的應(yīng)用前景，比如日本的下一代超級(jí)計(jì)算機(jī)富岳[13,14]就采用了SVE。

SVE引入了一些新的架構(gòu)特點(diǎn)，包括：(1)可變向量長(zhǎng)度：SVE提供32個(gè)向量寄存器，向量寄存器長(zhǎng)度為128的整數(shù)倍，最低128位，最高可支持2 048位；(2)每通道預(yù)測(cè)(Per-Lane Predication):SVE提供16個(gè)預(yù)測(cè)寄存器(Predicate Register)，預(yù)測(cè)寄存器每一位控制向量寄存器的一個(gè)字節(jié)，通過(guò)預(yù)測(cè)寄存器中每一位的有效狀態(tài)來(lái)控制向量寄存器中對(duì)應(yīng)元素是否參與運(yùn)算；(3)聚集取和分散存(Gather-Load，Scatter Store):支持非連續(xù)存儲(chǔ)數(shù)據(jù)的高效訪問(wèn)。

SVE提供了比傳統(tǒng)SIMD架構(gòu)更長(zhǎng)的向量寄存器，具備更強(qiáng)的計(jì)算能力，在SVE上運(yùn)行的程序理論上可獲得更高的性能。本文基于SVE提出了一種SPH的SIMD向量?jī)?yōu)化方法，利用SVE對(duì)SPH方法的關(guān)鍵部分粒子間相互作用的計(jì)算開(kāi)展向量化并行，探索SPH的數(shù)據(jù)級(jí)并行。實(shí)驗(yàn)結(jié)果表明，采用SVE可以獲得粒子間相互作用計(jì)算的顯著性能提升。

2 背景介紹

SVE是ARM在AArch64架構(gòu)基礎(chǔ)上針對(duì)高性能應(yīng)用所推出的SIMD擴(kuò)展指令集。ARM為SVE定義了32個(gè)向量寄存器和16個(gè)預(yù)測(cè)寄存器。向量寄存器的寬度是128的整數(shù)倍，范圍從128位～2 048位。預(yù)測(cè)寄存器的每一位(bit)對(duì)應(yīng)向量寄存器的一個(gè)字節(jié)(byte)，通過(guò)不同位的有效性控制向量寄存器中不同元素是否參與向量操作，預(yù)測(cè)寄存器可以看作是一個(gè)控制向量運(yùn)算的掩碼。ARM為SVE引入了VLA編程方式，對(duì)于循環(huán)的向量化可以實(shí)現(xiàn)與處理器向量寄存器寬度無(wú)關(guān)的控制。

2.1 每通道預(yù)測(cè)

圖1通過(guò)一條向量化融合乘加指令(fmla)簡(jiǎn)要說(shuō)明基于預(yù)測(cè)寄存器的向量操作。指令fmla的操作是把向量寄存器z1與z2中的數(shù)據(jù)相乘，乘積與向量寄存器z0中的數(shù)據(jù)相加，最后將結(jié)果寫(xiě)入向量寄存器z0中。向量寄存器的后綴·d表示寄存器中存放的是雙精度浮點(diǎn)(double)類(lèi)型的操作數(shù)。預(yù)測(cè)寄存器p0/m為合并模式，表示若預(yù)測(cè)寄存器該有效位無(wú)效，則目標(biāo)寄存器對(duì)應(yīng)通道數(shù)據(jù)保持原值。3個(gè)向量寄存器z0,z1和z2 中參與運(yùn)算的元素由預(yù)測(cè)寄存器p0對(duì)應(yīng)的位控制，若對(duì)應(yīng)的位為T(mén)RUE(1)，則向量寄存器中對(duì)應(yīng)的元素參與運(yùn)算，如A[7]、A[5]、A[3]、A[1]，反之若預(yù)測(cè)寄存器中對(duì)應(yīng)的位為FALSE(0)，則向量寄存器中對(duì)應(yīng)的元素不參與運(yùn)算，如A[6]、A[4]、A[2]、A[0]。

Figure 1 SIMD operation based on predicate register

2.2 VLA

以DAXPY運(yùn)算為例簡(jiǎn)要說(shuō)明VLA編程方式?；赩LA的DAXPY運(yùn)算代碼如下所示：

1 //-------------------------------------------------------------

2 // voiddaxpy(int64_tn,doubleda,double *dx,double *dy) {

3 // for (int64_ti=0;i

4 //dy[i]=dx[i] *da+dy[i];

5 // }

6 // }

7 //-------------------------------------------------------------

8 //x0=n,x1=&dy[0],x2=&dx[0],d0=da

9 daxpy:

10 movx3,0//x3=i=0

11 movz0.d,d0//z0=broadcast(da)

12whileltp0.d,xzr,x0//p0=while(0

13 ·L5:

14 ld1dz1.d,p0/z,[x1,x3,lsl 3]//p0?:z1=dy[i]

15 ld1dz2.d,p0/z,[x2,x3,lsl 3]//p0?:z2=dx[i]

16 fmadz1.d,p0/m,z0.d,z2.d/*p0?:z1 +=dx[i]*da*/

17 st1dz1.d,p0,[x2,x3,lsl 3]//p0?:dy[i]=z1

18 incdx3//i+=(VL/64)

19whileltp0.d,x3,x0//p0=while(i++

20 b.any.L5// more to do?

21ret

代碼中第2～6行顯示的是DAXPY運(yùn)算的源代碼，第9～21行顯示的是經(jīng)過(guò)編譯之后生成的含有SVE指令的匯編代碼。記處理器的向量寄存器寬度為VL，1次向量操作最多可處理的數(shù)據(jù)個(gè)數(shù)為VL/64。第10行和第11行初始化寄存器x3和z0。第12行的“whilelt”指令對(duì)預(yù)測(cè)寄存器p0初始化，分別對(duì)比0，1,…,j,…,VL/64-1和n的大小并為p0的對(duì)應(yīng)位賦值，即p0[j]=j

3 SVE向量化優(yōu)化方法

SPH將計(jì)算域內(nèi)的連續(xù)介質(zhì)離散成一系列具有質(zhì)量、速度和能量的“粒子”，粒子的物理量通過(guò)對(duì)其支持域內(nèi)其他粒子進(jìn)行插值近似，并將函數(shù)導(dǎo)數(shù)轉(zhuǎn)化為一個(gè)核函數(shù)的導(dǎo)數(shù)，從而實(shí)現(xiàn)對(duì)Navier-Stokes方程組的離散。粒子隨著流體流動(dòng)并與鄰近的粒子產(chǎn)生交互。SPH方法的計(jì)算過(guò)程被分解成一系列時(shí)間步，每一個(gè)時(shí)間步更新粒子的物理信息，其中最重要的計(jì)算過(guò)程是粒子間相互作用的計(jì)算。粒子間相互作用計(jì)算的串行算法如算法1所示。

算法1SPH粒子間相互作用計(jì)算串行算法

Procedure ParticlesInteractionSerial

Begin

1 // 記P為所有粒子的集合

2foreachp∈Pdo

3q←NeighborCellParticles(p);

4for(i=0;i

5d=distance(p,q[i]);

6if(d

7interact(p,q[i]);

8endfor

9endforeach

end

對(duì)于計(jì)算域P中的任意一個(gè)粒子p，首先通過(guò)函數(shù)NeighborCellParticles獲取該粒子鄰近單元中所有粒子的列表q，對(duì)于列表q中的任意一個(gè)粒子q[j]，根據(jù)p與q[j]的距離d判斷q[j]是否在粒子p的支持域內(nèi)部。若粒子q[j]在粒子p的支持域內(nèi)(d

函數(shù)interact通過(guò)如式(1)～式(3)所示的動(dòng)量方程、連續(xù)性方程和狀態(tài)方程計(jì)算粒子間相互作用。

動(dòng)量方程、連續(xù)性方程和狀態(tài)方程如式(1)～式(3)所示：

(1)

(2)

(3)

3.1 SPH粒子間相互作用向量化并行計(jì)算框架

在粒子間相互作用計(jì)算的串行算法中，內(nèi)層循環(huán)的每一次迭代都需要判斷2個(gè)粒子之間的距離是否小于支持域的半徑，即循環(huán)內(nèi)存在條件分支。條件分支是傳統(tǒng)SIMD優(yōu)化方法的一個(gè)難點(diǎn)，而SVE所提供的預(yù)測(cè)寄存器則能夠以類(lèi)似掩碼的形式很好地處理?xiàng)l件分支。

算法2描述了SPH粒子間相互作用計(jì)算的SVE向量化算法框架。

算法2SPH粒子間相互作用計(jì)算SVE向量化算法

Procedure ParticlesInteractionSVE

Begin

1 // 記P為所有粒子的集合

2foreachp∈Pdo

3q←NeighborCellParticles(p);

4for(i=0;i

5pd0=svwhilelt(i,len(q));

6vq[0:n-1]=ld(pd0,q[i]);

7vd[0:n-1]=distance(pd0,p,vq[0:n-1]);

8pd1=svcmplt(pd0,vd,h);

9interactSVE(pd1,p,vq[0:n-1]);

10endfor

11endforeach

end

如第4行所示，循環(huán)的每一次迭代最多處理n個(gè)粒子，n的計(jì)算方法為向量寄存器的寬度除以單個(gè)數(shù)據(jù)的大小(假如向量寄存器寬度為128位，數(shù)據(jù)類(lèi)型為雙精度浮點(diǎn)，即一個(gè)數(shù)據(jù)的大小為64位，則n=128/64=2)。第5行通過(guò)循環(huán)控制變量j和列表q的長(zhǎng)度len(q)構(gòu)造預(yù)測(cè)寄存器pd0，預(yù)測(cè)寄存器pd0包含n個(gè)元素，pd0[k]=(j+k)

觀察粒子間相互作用計(jì)算的3個(gè)方程，動(dòng)量方程和連續(xù)性方程均涉及大量的求和運(yùn)算，占據(jù)了粒子間相互作用的絕大部分計(jì)算時(shí)間，而且2個(gè)方程中存在大量的數(shù)據(jù)級(jí)并行，因此，下面將重點(diǎn)描述如何采用SVE對(duì)動(dòng)量方程和連續(xù)性方程進(jìn)行向量化。以下的向量化并行計(jì)算均在預(yù)測(cè)寄存器的控制下進(jìn)行。

3.2 動(dòng)量方程向量化

動(dòng)量方程是關(guān)于粒子速度矢量vk的方程，包含x,y,z3個(gè)方向的分量，各對(duì)應(yīng)一個(gè)方程，如式(4)～式(6)所示：

(4)

(5)

(6)

Figure 2 SIMD calculation of

動(dòng)量方程右側(cè)的第2部分比第1部分更為復(fù)雜，將分3步描述第2部分的向量化計(jì)算方法。首先在圖3中描述了rij·vij的計(jì)算過(guò)程，其中，rij·vij=rij|x×vij|x+rij|y×vij|y+rij|z×vij|z,圖3中記rvij=rij·vij。注意圖3中(rvij0|x,rvij1|y,…,rvijn-1|y),(rvij0|z,rvij1|z,…,rvijn-1|z)與(rvij0|x,rvij1|x,…,rvijn-1|x)的計(jì)算過(guò)程類(lèi)似，故省去。在圖3中，ri|x和vi|x通過(guò)標(biāo)量取獲得，(rj0|x,rj1|x,…,rjn-1|x)和(vj0|x,vj1|x,…,vjn-1|x)通過(guò)向量取獲得，圖3中的計(jì)算全部為向量運(yùn)算。

Figure 3 SIMD calculation of rij·vij

Figure 4 SIMD calculation of

記:

Figure 5 SIMD calculation of

3.3 連續(xù)性方程向量化

Figure 6 SIMD calculation of wij

4 實(shí)驗(yàn)結(jié)果

本節(jié)對(duì)本文所提出的向量化方法的加速效果進(jìn)行評(píng)估。本文采用5個(gè)測(cè)試用例，分單精度和雙精度2種模式，對(duì)比了從128位～2 048位共16種向量寄存器寬度情況下粒子間相互作用向量化算法對(duì)串行算法的加速比，得到了單精度和雙精度2種模式下最高3.71倍和4.92倍的加速比。

(1)模擬器。本文采用了ARM指令模擬器ARMIE(ARM Instruction Emulator)對(duì)SVE指令開(kāi)展模擬。ARMIE是ARM公司為了方便相關(guān)從業(yè)人員進(jìn)行SVE指令的研究而推出的一個(gè)模擬器，ARMIE可以通過(guò)選項(xiàng)“-msve-vector-bits”指定模擬的向量寄存器寬度，支持從128位～2 048位共16種不同的向量寄存器寬度。受模擬器運(yùn)行速度的限制，對(duì)測(cè)試用例完整計(jì)算過(guò)程的模擬難以在有限時(shí)間內(nèi)完成，因此后續(xù)實(shí)驗(yàn)針對(duì)每一個(gè)測(cè)試用例選取了計(jì)算過(guò)程中的部分時(shí)間迭代開(kāi)展模擬。

(2)軟硬件平臺(tái)。ARMIE需要運(yùn)行在支持ARM AArch64指令集的硬件平臺(tái)上。實(shí)驗(yàn)所采用的是一個(gè)FT-2000+/64[15]的CPU，包含64個(gè)核，64 GB主存，實(shí)驗(yàn)使用的編譯器為GCC-10.1.0，為了生成SVE指令，給編譯器提供的編譯選項(xiàng)為“-O2-march=armv8-a+sve”。

(3)SPH軟件平臺(tái)。本文選擇開(kāi)源SPH項(xiàng)目DualSPHysics[16]作為本實(shí)驗(yàn)的軟件平臺(tái)。DualSPHysics由西班牙維戈大學(xué)和英國(guó)曼徹斯特大學(xué)共同開(kāi)發(fā)，采用C++語(yǔ)言編寫(xiě)，適用于劇烈運(yùn)動(dòng)或重力場(chǎng)和液體的建模場(chǎng)景。本文對(duì)DualSPHysics項(xiàng)目中計(jì)算粒子間相互作用的部分進(jìn)行了向量?jī)?yōu)化，粒子間相互作用的計(jì)算占據(jù)了DualSPHysics計(jì)算過(guò)程的大部分時(shí)間。對(duì)于粒子間相互作用的加速可以顯著提高整個(gè)軟件的性能。

(4)測(cè)試用例。實(shí)驗(yàn)采用5個(gè)測(cè)試用例來(lái)驗(yàn)證向量化的加速效果，對(duì)測(cè)試用例的簡(jiǎn)要描述以及測(cè)試用例的規(guī)模(粒子數(shù)量)如表1所示。

對(duì)于上述5個(gè)測(cè)試用例，本文采用了2種不同的精度模式：?jiǎn)尉群碗p精度，對(duì)比了從128位～2 048位16種不同向量寄存器寬度的情況下，粒子間相互作用計(jì)算部分，經(jīng)過(guò)SVE向量化之后對(duì)比非向量化的加速比。

Table 1 Description of the testcases

單精度和雙精度2種不同模式下的加速比如圖7和圖8所示。可以看出雙精度模式下的加速比要高于單精度模式下的加速比。單精度模式和雙精度模式下的最高加速比來(lái)自于測(cè)試用例Pump，分別為3.71倍和4.92倍。加速效果與ARM針對(duì)其他Benchmark所給出的最高加速比存在差距的原因主要包括2個(gè)方面：一方面是因?yàn)榱Ｗ娱g相互作用的計(jì)算仍然存在串行計(jì)算部分，并不是所有的計(jì)算都可以通過(guò)向量化并行計(jì)算，比如由于粒子的位置發(fā)生了變化，需要更新粒子所在單元的信息；另外一方面，在算法2的SVE向量化算法中，第4行所遍歷的粒子均需在第7行計(jì)算其與粒子的距離，只有距離小于支持域半徑的粒子才會(huì)參與后續(xù)的向量運(yùn)算，向量寄存器并沒(méi)有被填滿，這也是SVE的特點(diǎn)之一。

Figure 7 Speedup in single precision mode

Figure 8 Speedup in double precision mode

為了驗(yàn)證SVE向量化算法的正確性，本文以測(cè)試用例DamBreakCubes為例，對(duì)比了非向量化算法的渲染結(jié)果與SVE向量化算法的粒子分布，結(jié)果如圖9和圖10所示。通過(guò)對(duì)比可以發(fā)現(xiàn)，經(jīng)過(guò)SVE向量化加速之后的粒子分布與原始結(jié)果的渲染圖一致，表明本文所提出的SVE向量化方法是正確的。

Figure 9 Rendering effect of the testcase DamBreakCubes using serial algorithm

Figure 10 Particle distribution of the testcase CamBreakCubes using SVE algorithm

5 結(jié)束語(yǔ)

本文基于ARM SVE擴(kuò)展指令集對(duì)光滑粒子流體動(dòng)力學(xué)(SPH)開(kāi)展了SIMD優(yōu)化，對(duì)SPH計(jì)算過(guò)程中的粒子間相互作用計(jì)算進(jìn)行了向量化并行，在單精度模式下取得了最高3.71倍的加速比，雙精度模式下取得了4.92倍的加速比。下一步將開(kāi)展大規(guī)模SPH計(jì)算在SVE硬件平臺(tái)上的加速效果研究。