基于ARM SVE的光滑粒子流體動力學(xué)SIMD加速方法*

2021-06-25 09:46:04范小康夏澤宇龍思凡楊燦群

計算機工程與科學(xué) 2021年6期

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

(國防科技大學(xué)計算機學(xué)院,湖南長沙 410073)

1 引言

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

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

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

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

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

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

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

2 背景介紹

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

2.1 每通道預(yù)測

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

Figure 1 SIMD operation based on predicate register

2.2 VLA

以DAXPY運算為例簡要說明VLA編程方式?；赩LA的DAXPY運算代碼如下所示：

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

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

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

算法1SPH粒子間相互作用計算串行算法

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

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

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

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

(1)

(2)

(3)

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

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

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

算法2SPH粒子間相互作用計算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個粒子，n的計算方法為向量寄存器的寬度除以單個數(shù)據(jù)的大小(假如向量寄存器寬度為128位，數(shù)據(jù)類型為雙精度浮點，即一個數(shù)據(jù)的大小為64位，則n=128/64=2)。第5行通過循環(huán)控制變量j和列表q的長度len(q)構(gòu)造預(yù)測寄存器pd0，預(yù)測寄存器pd0包含n個元素，pd0[k]=(j+k)

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

3.2 動量方程向量化

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

(4)

(5)

(6)

Figure 2 SIMD calculation of

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 實驗結(jié)果

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

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

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

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

(4)測試用例。實驗采用5個測試用例來驗證向量化的加速效果，對測試用例的簡要描述以及測試用例的規(guī)模(粒子數(shù)量)如表1所示。

對于上述5個測試用例，本文采用了2種不同的精度模式：單精度和雙精度，對比了從128位～2 048位16種不同向量寄存器寬度的情況下，粒子間相互作用計算部分，經(jīng)過SVE向量化之后對比非向量化的加速比。

Table 1 Description of the testcases

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

Figure 7 Speedup in single precision mode

Figure 8 Speedup in double precision mode

為了驗證SVE向量化算法的正確性，本文以測試用例DamBreakCubes為例，對比了非向量化算法的渲染結(jié)果與SVE向量化算法的粒子分布，結(jié)果如圖9和圖10所示。通過對比可以發(fā)現(xiàn)，經(jīng)過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é)束語

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