激光等離子體相互作用模擬的并行和加速研究*

武海鵬，文敏華+，SEE Simon，林新華,3

1.上海交通大學 高性能計算中心，上海 200240

2.NVIDIA Technology Center，新加坡

3.東京工業大學 學術國際情報中心，日本 東京

1　問題簡介

近20年來，生成超短激光脈沖能力的不斷增強極大地促進了激光等離子體相互作用領域的理論和實驗研究[1-2]。激光驅動的電子加速是其中最重要的研究方向之一[3]。

PIC（particle-in-cell）方法已經被廣泛應用在激光等離子體相互作用的模擬和其他物理學模擬中[4]。這種模擬方法結果的質量依賴于系統中引入的大量粒子，因此基于PIC方法的應用會包含大量的計算。由于可能存在大量的數據沖突和不規則內存訪問，如何很好地將應用并行化并且利用好集群的特性來加速程序是一個很有挑戰性的任務。

本文研究用到的代碼“2DPIC”是一個新研發的電磁學的基于PIC方法的激光等離子體相互作用模擬的代碼[5]。在這個代碼中，應用了方向劃分方法來求解麥克斯韋方程組。不同于有限差分時間域的方法，這個方法不受限于Courant條件：Δt=Δx/c<Δy/c，其中Δt、Δx和Δy分別是時間、縱向和橫向分辨率[6]。

在將程序移植到GPU上的時候，利用“數組”作為存儲網格和粒子數據的主要數據結構。為了處理“Scatter”階段中存在的數據沖突，在最開始嘗試利用了傳統的原子操作的方法。然而發現大量的原子操作成為程序性能的瓶頸，因此提出了動態冗余算法和混合精度算法來加速這部分計算過程。

另外值得注意的一點是，當程序運行在安裝了很多GPU的集群中時，MPI（message passing interface）通信部分的運行時間會顯著增加。在這種情況下，利用 GPUDirect RDMA（remote direct memory access）技術可以將數據傳輸過程加速30%～60%，取決于要傳輸到其他MPI區域的數據量的多少。

2　相關工作

近些年，盡管一些基于PIC方法的代碼已經被移植到了GPU平臺，但是本文進行了一些不同的探索。

Decyk和Singh將一個比較簡單的PIC代碼移植到了GPU，這個代碼只應用了Poisson方程來計算磁場信息[7]，但是本文用到的代碼是基于電磁學，并且有更大的計算強度。另一方面，本文所用的粒子排序算法不同，在實現中，粒子在每一步迭代后都會保持有序，而本文參數化了這個過程。對于實驗方面，他們利用了一塊Tesla C1060和GTX 280，本文在一個GPU集群中進行了實驗評估。

Joseph等人將PIC算法中計算強度最大的部分在一個GPU上進行了并行化[8]，但是本文將整個算法都移植到了GPU上。他們描述了一種新的并行的分層的木桶排序算法，其中粒子在每次循環迭代后逐漸變得有序。

Chen等人描述了一種基于CPU+GPU的一維PIC算法的移植實現[9]。這種方法的基本特點是，從場求解器中分離出了粒子軌跡的計算部分，于此同時也保持了系統自洽。

Kong等人研究了一種線程競爭的算法來解決數據沖突的問題，但是這種算法仍然類似于傳統的單精度的原子操作方法，并且他們只是在一塊GPU上進行了測試和分析[10]。

Burau等人研究了一種可擴展到集群中的針對等離子體物理模擬的PIC方法的實現[11]。他們利用了LinkedList作為主要的數據結構，并且在4個GPU上進行了測試。對于粒子排序，他們通過實驗發現了在進行一定次數的循環之后對粒子進行重排序可以加速程序，但沒有詳細分析并且給出一個解決方案。

Wang等人開發了基于PIC方法的GTC-P代碼，并且在一些超級計算機上進行了測試[12]。這個代碼的算法和激光等離子體相互作用的模擬算法是不同的，他們所用的一些優化手段不適合應用在本文的模擬代碼中。另一方面，他們采用傳統的MPISend-Recv模式在CPU和GPU之間交換數據，沒有利用GPUDirect RDMA技術。

因此和上述方法相比，本文采用的實現和優化方法有以下的不同：

（1）在整體上把模擬算法移植到了GPU上，而不僅僅是其中計算強度高的部分。

（2）探索了一系列方法來減少由于原子操作而產生的內存沖突帶來的性能下降。同時也考慮到了盡量減少內存的使用這一因素。

（3）研究了一種參數化的粒子排序算法，不需要每次迭代之后都進行粒子排序，從而在取得較好性能的同時不會浪費過多的GPU內存。

（4）探索了在激光等離子體相互作用模擬中應用GPUDirect RDMA技術，而不是利用傳統的MPISendRecv方法在CPU和GPU之間進行數據傳輸。

3　算法

這個代碼同時完成對場信息求解的麥克斯韋方程組和對微粒子的運動求解的方程組的計算。對于沒有發生碰撞的等離子體，在每一個時間步Δt內，下列相對論方程組會作用于每一個粒子：

粒子會對電荷和電流密度生成貢獻值。在當前網格上會求解以下麥克斯韋方程組：

這個過程會不斷迭代進行，直到在等離子體場中達到自洽的狀態。PIC循環如圖1所示，主要包含4個階段。

Fig.1　PIC cycle圖1　PIC循環邏輯

（1）Gather階段：粒子的電場和磁場信息是由粒子在網格中的相對位置決定的網格點處得到的。每個粒子通過這些貢獻值來生成當前位置的力。

（2）Push階段：上一個Gather階段產生的力推動粒子到新的位置。

（3）Scatter階段：每一個粒子將自己的貢獻傳播到當前網格中。每一個網格點上的貢獻值被累積起來加到局部密度中。

（4）Field階段：每一個網格點會從相鄰的網格點來取得數據計算出新的電場和磁場值。

4　基于GPU的并行化

4.1　合并函數

每一個PIC方法的循環都由圖1描述的幾個主要階段組成。盡管每次迭代可以看作依賴于前一次迭代的結果（粒子新的位置），每一個粒子的gather、push和scatter這3個階段可以看作是相互獨立的，因此把這3個階段合并為一個kernel，叫作“ParticleKernel”。通過這種方法，不僅簡化了代碼，而且更為重要的是利用了數據的局部性，這意味著對于每一個粒子，在上一個階段產生的數據將馬上被用在下一個階段的計算中。偽代碼如下所示：

4.2　線程分配策略

為全局粒子數組中的每一個位置都分配了一個GPU線程。首先判斷對應于這個粒子的位置是否“有效”。然后利用上一個階段的輸出作為下一個階段的輸入依次完成Gather、Push和Scatter階段。

4.3　數據結構轉換

PIC方法主要基于兩種類型的數據集合：一種是坐標連續的充電粒子相關的數據，另一種是網格相關的數據。對于每一種數據，都額外分配一個全局數組來存儲。通過這種從AOS（array-of-structure）到SOA（structure-of-array）的轉換，可以最大限度地減少數據量，無論是通過CudaMemcpy從CPU傳到GPU的數據量，還是通過MPI通信在不同區域中傳輸的數據量。

5　對GPU版本的優化

5.1　“Scatter”階段優化

如下面偽代碼所示，“ParticleKerel”一共有3個步驟：計算貢獻值，確定新的目標地址，累加貢獻值。

如圖2所示，每一個粒子都可能向9個方向移動（8個鄰接位置和當前位置），并且可能會對16個位置（圖中的灰色區域）進行寫操作。因此第三步“將貢獻值加到當前網格”會帶來大量的數據沖突，從而將嚴重地降低程序性能。為了解決這個問題，最簡單的方法是利用傳統的“原子操作”方法。然而由于CUDA現在并沒有提供基于雙精度的原子操作的底層實現，本文嘗試應用CAS（compare&set）方法來實現。但是經過分析程序，發現這個部分仍然占用了超過80%的時間，因此首先分析原因，然后依次介紹本文的優化方法。

Fig.2　Data conflicts on grid圖2　粒子產生數據沖突

5.1.1　動態冗余算法

在嘗試分析了產生數據沖突的根本原因之后，總結了兩方面因素：一是同屬一個網格的粒子可能嘗試去寫同一個內存地址；另一方面，屬于不同網格的粒子也可能產生數據寫沖突。比較直接的想法是可以分配一塊和原始的網格大小相同的內存區域來存儲由于數據沖突而產生的臨時數據。

更進一步，本文決定采用兩個層次的“數據冗余”方法。一方面，如圖3所示，為那些由于向不同的網格嘗試寫操作而產生數據沖突的粒子分配了冗余空間。因為一個粒子可能向16個方向進行寫操作，所以對整個網格復制了16份。接下來考慮每一個粒子，原來對“原始”數組的寫操作現在會轉移到另外的冗余數組中，取決于數據會寫入16個方向中的哪一個網格。通過這種方法，原來由于寫入不同的網格而產生的數據寫操作將會被分割開，從而大大減少了這種情況下的數據沖突數量。

Fig.3　Data conflicts for particles in different cells圖3　解決屬于不同cell的粒子產生的數據沖突

另一方面，那些屬于同一個網格中的粒子也可能產生數據沖突，如圖4所示，對這一個網格進行了冗余化，將當前網格的粒子按照它們的位置進行了劃分。具體來說，首先將一個網格劃分為多份，然后對每一個粒子計算并判斷它們屬于哪一個子區域。之后對當前網格信息的更新操作會在冗余數組上進行。通過這種方法，分割了屬于同一個網格可能產生的數據寫操作。為了控制在這種情況下為每一個網格分配多少冗余空間，引入了一個參數duplication_num。通過仔細地選擇這個參數合適的值，可以找到在減少數據沖突和管理數據的開銷之間的平衡點。

Fig.4　Data conflicts for particles in same cells圖4　解決屬于同一個cell的粒子產生的數據沖突

像上文提到的，當所有的粒子寫操作完成時，會從所有的冗余數組中收集這些臨時數據，并在另一個Kernel函數中將值更新到原始數組中。然而這種方法很有可能會占用太多的內存資源。為了解決這個問題，本文引入了“動態冗余”的思想，偽代碼如下所示。這個算法的基本思想是考慮到粒子在整個網格上其實不是均勻分布的，可以把粒子的分布情況考慮在內。不是簡單地為每一個網格都分配相同數目的冗余空間，而是分配給那些包含粒子數目多的網格更多的冗余空間，相反，給那些包含粒子較少的網格較少的冗余空間。具體視硬件資源的情況而定。這個方法可以極大地減少內存空間的使用，尤其是當粒子分布很不均勻的情況。當程序運行在28個GPU的集群中時，與初始的GPU版本相比，取得了約1.7倍的加速比。

現在可以從整體上來分析“動態冗余算法”。從圖5可以看出，從兩個維度減少了數據沖突：一個是從各個方向生成冗余空間，另一個是為每一個網格生成復制。“動態”意味著每一個子網格中的冗余空間的數量是不同的，取決于粒子在整個網格上分布的“密度”。

5.1.2　混合精度算法

另外一個來優化“Scatter”階段的方法是加速原子操作。盡管已經用傳統的CAS方法來實現雙精度的原子操作，但是這種方法仍然會比CUDA運行時環境提供的單精度原子方法慢得多。為了能擁有單精度運算的快速和雙精度運算的精確，本文將二者結合了起來。

具體來說就是在ParticleKernel之前和之后利用另外兩個GPU Kernel來做單精度和雙精度之間的轉換。代碼的其余部分，如數據的初始化，Field階段和診斷階段等，仍然用雙精度來進行計算，盡量減少由于精度損失帶來的影響。本文應用這個方法在結果精度可以接受的情況下取得了約2倍的加速比。

Fig.5　Duplications in two dimensions圖5　在兩個維度上通過冗余解決數據沖突

5.2　粒子排序算法

粒子可能逃逸出整個網格區域或者移動到相鄰的MPI域中。在這兩種情況中，全局數組中存儲這個粒子的位置都會變成“空”或者“無效”的位置。如果不去處理，那么這些位置就會變得越來越多。這樣不僅會浪費內存資源，也會浪費GPU計算資源，因為這里采用的GPU線程分配策略是為每一個全局粒子數組的位置分配一個GPU線程。進行粒子排序的另一個目的是使得在物理位置上鄰近的粒子在內存中也能連續存儲。盡管排序粒子能帶來較好的性能，但是同時排序也會耗費一定的時間和計算資源，在有些情況下，過于頻繁地進行粒子排序并不能有效提升程序性能。

為了在排序帶來的好處和產生的額外開銷之間取得平衡，引入另外一個可調節的參數empty_ratio，這個參數表示了當前粒子數組的“空余程度”。如果數組中空余位置數量和數組的總位置數量的比值大于本文設置的empty_ratio，那么粒子就會根據它們在整個網格中的位置進行排序，如圖6所示。通過調節這個參數值，可以回收“無效”的空間，從而在得到較好的程序性能的同時，也能最大程度地節省內存空間。

Fig.6　Particle reordering圖6　粒子排序

5.3　CUDA-Aware MPI和GPUDirect RDMA

MPI能很好地和CUDA編程模型相兼容[13]。對于一般的MPI實現來說，只有指向host空間的指針可以被當作參數傳遞給MPI函數。然而如果把MPI和CUDA相結合，就可以發送GPU的緩存，而不僅僅是CPU端的數據了。如果不利用CUDA-Aware的MPI實現，就需要利用cudaMemcpy函數把GPU內存中的數據傳輸到CPU端的內存。然而如果利用了這個技術，GPU端的緩存可以直接被MPI傳輸到另一端。從下面的代碼對比中可以看出，使用了CUDA-Aware的MPI實現后，代碼會變得更簡潔，編程也更加容易。

除了更高的可用性之外，CUDA-Aware MPI還有哪些優點呢？它不僅使得MPI+CUDA更容易使用，而且也使得應用程序取得更好的運行效率，這基于以下兩點原因[14]：

（1）所有進行數據傳輸的操作都會被流水化。

（2）MPI庫可以透明化地使用一些如GPUDirect的加速技術。

遠程DMA（RDMA）技術是在CUDA5.0中引入的GPUDirect技術中重要的一部分[15]，它連通了GPU和第三方應用了PCI-E標準的硬件。GPU緩沖區中的數據可以不通過CPU端而直接被送到網卡進行傳輸，從而消除了CPU到GPU和其他PCI-E設備的內存帶寬占用。這樣也就顯著增大了GPU和其他節點的MPISendRecv的效率。在一個由許多GPU加速節點組成的集群中，嘗試應用了這個技術來加速數據傳輸過程。通過實驗發現，當傳輸的數據量比較大時會有比較好的加速效果。隨著將來越來越多的加速卡和節點的投入使用，更好地利用這項技術來減少數據傳輸的消耗將會變得越來越重要。

6　實驗結果和分析

表1展示了本文的實驗環境。在這個模擬中，設置 Δx=λ/200，Δy=λ/100和 Δt=Δx/c，在每一個非真空的網格中，對于電子、氘核和氚核，初始分別分配了64、32和32個粒子。這個模擬從t=-λ/C開始，與此同時激光脈沖從左邊界x=-λ處發射。

Table 1　Software and hardware environment表1　實驗的軟硬件環境

從圖7可以看出，總體來說，如果把使用單個GPU的初始CUDA版本作為對比，在使用這些優化方法后，與初始版本相比取得了更好的可擴展性。下面分析每一種優化方法。

Fig.7　Speedup after utilizing all optimizations圖7　利用所有優化方法之后的加速比

6.1　動態冗余算法

本文測試了利用動態冗余算法之后與初始的GPU版本使用相同數量GPU相比的加速情況。從圖8可以看出，當用兩個GPU時加速比是2.18。隨著使用的GPU數量的增多，加速比會緩慢下降。當使用28塊GPU時，加速比為1.69。

Fig.8　Speedup when utilizing different optimizations(with same number of GPUs)圖8　使用各種優化方法之后的加速比（使用相同數量GPU)

這里不考慮MPI通信時間隨著GPU數量增多而變化帶來的影響，因為使用冗余優化方法并不會改變通過MPI傳輸的數據數量。由于總的粒子數目是固定的，隨著使用GPU數的增加，每個GPU上的粒子數目會減少，這樣數據沖突的數量就會降低，冗余優化的效果不突出，反而從冗余空間收集數據和更新數組的額外開銷變得明顯，導致加速比緩慢下降。

6.2　混合精度算法

總體上來說，隨著GPU數目的增多，混合精度計算方法帶來的加速收益會變得越來越明顯，當GPU數量超過8后，加速效果會超過動態冗余算法。這是因為雖然總的計算量是固定的，但是當GPU數目較少時，每個GPU上的計算量就會較多，這樣進行數據精度轉換帶來的額外開銷會比較大。當每個GPU上的粒子數變少后，由于各個GPU上精度轉換函數是并行執行的，精度轉換的總開銷也會相應降低，加速效果就會變得越來越明顯。

6.3　RDMA

對于GPUDirect RDMA，本文通過改變MPI線程的數目來改變需要傳輸的數據量。發現使用GPUDirect RDMA可以顯著減少數據傳輸的時間。如圖9所示，當傳輸的數據量大于3 KB時，可以獲得2.8倍的加速比。但是當傳輸的數據量為0.9 KB時，加速比只有1.37。這是因為當需要傳輸的數據量比較少時，RDMA的優勢不明顯。待傳輸的數據量越大，越多的數據要傳輸到另一個MPI域中。對于傳統的MPISendRecv方法來說，就需要更多的時間將數據從CPU端傳輸到GPU端。

Fig.9　Speedup for MPI communication part圖9　MPI通信時間的加速比

這里需要指出的一點是，對于非RDMA版本，MPI通信時間的計算不僅僅包括純粹的數據傳輸時間，還包括收集數據，將數據放到數組中和把接收到的數據更新到正確的數組位置中的時間，因為RDMA版本并不需要這些操作。

6.4　粒子排序

本文的粒子排序算法的中心思想是引入了一個參數empty_ratio，目的是在性能達到或者接近最優的情況下找到這個參數的最小值。用8個MPI進程來測試這個參數不同取值情況下的程序性能，結果如圖10所示。實驗發現當empty_ratio小于0.05時，性能受這個參數的影響很大。然而當大于0.05時，性能接近于穩定。因為這個參數越小代表著內存空間的使用越少，所以可以通過設置這個參數為0.05來達到性能和內存使用之間的平衡。

Fig.10　Performance when setting different empty_ratio圖10　不同的空余率對程序性能的影響

7　結論

本文研究了一種把整體的基于PIC方法的激光等離子體相互作用模擬的代碼移植到GPU端的可行方法，并且對比原始的CPU代碼取得了可觀的加速比。基于這個初始的GPU版本，本文介紹了一系列的優化方法來加速基于原子操作的包含大量數據沖突的Scatter階段，包括動態冗余算法、混合精度計算方法和一種參數化的粒子排序方法。也嘗試利用并且評估了GPUDirect RDMA方法在集群中加速MPI通信時間的效果。發現當數據傳輸量大于一定閾值的情況下，這種方法能顯著減少MPI通信時間。相信這些優化方法也能應用于其他基于PIC方法的物理學模擬代碼，并且也對激光等離子體相互作用的研究發展有著非常重要的意義。

致謝林新華特別致謝日本學術振興會JSPS的RONPAKU項目資助。感謝NVIDIA GCOE的支持。感謝上海交通大學物理與天文系的翁蘇明教授的指導和幫助。

[1]Mourou G,Tajima T,Bulanov S.Optics in the relativistic regime[J].Reviews of Modern Physics,2006,78(2):309-371.

[2]Korzhimanov AV,Gonoskov AA,Khazanov E A,et al.Horizons of petawatt laser technology[J].Physics-Uspekhi,2011,54(1):9-28.

[3]Kostyukov I Y,Pukhov A M.Plasma-based methods for electron acceleration:current status and prospects[J].Physics-Uspekhi,2015,58(1):81-88.

[4]Wen Minhua,Yu Zhanpeng,See S,et al.A NVIDIA Kepler based acceleration of PIC method[J].Procedia Engineering,2013,61:398-401.

[5]Weng S M,Murakami M,Azechi H,et al.Quasi-monoenergetic ion generation by hole-boring radiation pressure acceleration in inhomogeneous plasmas using tailored laser pulses[J].Physics of Plasmas,2014,21:012705.

[6]Pfund R E W,Lichters R,Meyer-ter-Veh J.LPIC++a parallel one-dimensional relativistic electromagnetic particle-in-cell code for simulating laser-plasma-interaction[C]//Proceedings of the International Conference on Superstrong Fields in plasmas,Varenna,Aug 27-Sep 2,1997.

[7]Decyk V K,Singh T V.Adaptable particle-in-cell algorithms for graphical processing units[J].Computer Physics Communications,2011,182(3):641-648.

[8]Joseph R G,Ravunnikutty G,Ranka S,et al.Efficient GPU implementation for particle in cell algorithm[C]//Proceedings of the 25th International Symposium on Parallel and Distributed Processing,Anchorage,May 16-20,2011.Piscataway:IEEE,2011:395-406.

[9]Chen G,ChacóN L,Barnes D C.An efficient mixed-precision,hybrid CPU-GPU implementation of a non-linearly implicit one-dimensional particle-in-cell algorithm[J].Journal of Computational Physics,2012,231(16):5374-5388.

[10]Kong Xianglong,Huang M C,Ren Chuang,et al.Particle-incell simulations with charge-conserving current deposition on graphical processing units[J].Journal of Computational Physics,2011,230(4):1676-1685.

[11]Burau H,Widera R,Ho?Nig W,et al.PIConGPU:a fully relativistic particle-in-cell code for a GPU cluster[J].IEEE Transactions on Plasma Science,2010,38(10):2831-2839.

[12]Wang Bei,Ethier S,Tang W,et al.Modern Gyrokinetic particlein-cell simulation of fusion plasmas on top supercomputers[J].Computing Research Repository,arXiv:1510.05546,2015.[13]Kraus J.An introduction to CUDA-aware MPI[EB/OL].(2013-03-13).https://devblogs.nvidia.com/parallelforall/introductioncuda-aware-mpi/.

[14]Rossetti D.Benchmarking GPUDirect RDMA on modern server platforms[EB/OL].(2014-10-07).https://devblogs.nvidia.com/paralle-for-all/benchmarking-gpudirect-rdma-on-modernserver-platforms/.

[15]Corporation N.Developing a Linux kernel module using GPUDirect RDMA[EB/OL].(2017-06-23).http://docs.nvidia.com/cuda/gpudirect-rdma.