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

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

1.上海交通大學(xué) 高性能計(jì)算中心，上海 200240

2.NVIDIA Technology Center，新加坡

3.東京工業(yè)大學(xué) 學(xué)術(shù)國際情報(bào)中心，日本 東京

1　問題簡介

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

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

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

在將程序移植到GPU上的時(shí)候，利用“數(shù)組”作為存儲網(wǎng)格和粒子數(shù)據(jù)的主要數(shù)據(jù)結(jié)構(gòu)。為了處理“Scatter”階段中存在的數(shù)據(jù)沖突，在最開始嘗試?yán)昧藗鹘y(tǒng)的原子操作的方法。然而發(fā)現(xiàn)大量的原子操作成為程序性能的瓶頸，因此提出了動態(tài)冗余算法和混合精度算法來加速這部分計(jì)算過程。

另外值得注意的一點(diǎn)是，當(dāng)程序運(yùn)行在安裝了很多GPU的集群中時(shí)，MPI（message passing interface）通信部分的運(yùn)行時(shí)間會顯著增加。在這種情況下，利用 GPUDirect RDMA（remote direct memory access）技術(shù)可以將數(shù)據(jù)傳輸過程加速30%～60%，取決于要傳輸?shù)狡渌鸐PI區(qū)域的數(shù)據(jù)量的多少。

2　相關(guān)工作

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

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

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

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

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

Burau等人研究了一種可擴(kuò)展到集群中的針對等離子體物理模擬的PIC方法的實(shí)現(xiàn)[11]。他們利用了LinkedList作為主要的數(shù)據(jù)結(jié)構(gòu)，并且在4個(gè)GPU上進(jìn)行了測試。對于粒子排序，他們通過實(shí)驗(yàn)發(fā)現(xiàn)了在進(jìn)行一定次數(shù)的循環(huán)之后對粒子進(jìn)行重排序可以加速程序，但沒有詳細(xì)分析并且給出一個(gè)解決方案。

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

因此和上述方法相比，本文采用的實(shí)現(xiàn)和優(yōu)化方法有以下的不同：

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

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

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

（4）探索了在激光等離子體相互作用模擬中應(yīng)用GPUDirect RDMA技術(shù)，而不是利用傳統(tǒng)的MPISendRecv方法在CPU和GPU之間進(jìn)行數(shù)據(jù)傳輸。

3　算法

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

粒子會對電荷和電流密度生成貢獻(xiàn)值。在當(dāng)前網(wǎng)格上會求解以下麥克斯韋方程組：

這個(gè)過程會不斷迭代進(jìn)行，直到在等離子體場中達(dá)到自洽的狀態(tài)。PIC循環(huán)如圖1所示，主要包含4個(gè)階段。

Fig.1　PIC cycle圖1　PIC循環(huán)邏輯

（1）Gather階段：粒子的電場和磁場信息是由粒子在網(wǎng)格中的相對位置決定的網(wǎng)格點(diǎn)處得到的。每個(gè)粒子通過這些貢獻(xiàn)值來生成當(dāng)前位置的力。

（2）Push階段：上一個(gè)Gather階段產(chǎn)生的力推動粒子到新的位置。

（3）Scatter階段：每一個(gè)粒子將自己的貢獻(xiàn)傳播到當(dāng)前網(wǎng)格中。每一個(gè)網(wǎng)格點(diǎn)上的貢獻(xiàn)值被累積起來加到局部密度中。

（4）Field階段：每一個(gè)網(wǎng)格點(diǎn)會從相鄰的網(wǎng)格點(diǎn)來取得數(shù)據(jù)計(jì)算出新的電場和磁場值。

4　基于GPU的并行化

4.1　合并函數(shù)

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

4.2　線程分配策略

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

4.3　數(shù)據(jù)結(jié)構(gòu)轉(zhuǎn)換

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

5　對GPU版本的優(yōu)化

5.1　“Scatter”階段優(yōu)化

如下面?zhèn)未a所示，“ParticleKerel”一共有3個(gè)步驟：計(jì)算貢獻(xiàn)值，確定新的目標(biāo)地址，累加貢獻(xiàn)值。

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

Fig.2　Data conflicts on grid圖2　粒子產(chǎn)生數(shù)據(jù)沖突

5.1.1　動態(tài)冗余算法

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

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

Fig.3　Data conflicts for particles in different cells圖3　解決屬于不同cell的粒子產(chǎn)生的數(shù)據(jù)沖突

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

Fig.4　Data conflicts for particles in same cells圖4　解決屬于同一個(gè)cell的粒子產(chǎn)生的數(shù)據(jù)沖突

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

現(xiàn)在可以從整體上來分析“動態(tài)冗余算法”。從圖5可以看出，從兩個(gè)維度減少了數(shù)據(jù)沖突：一個(gè)是從各個(gè)方向生成冗余空間，另一個(gè)是為每一個(gè)網(wǎng)格生成復(fù)制?！皠討B(tài)”意味著每一個(gè)子網(wǎng)格中的冗余空間的數(shù)量是不同的，取決于粒子在整個(gè)網(wǎng)格上分布的“密度”。

5.1.2　混合精度算法

另外一個(gè)來優(yōu)化“Scatter”階段的方法是加速原子操作。盡管已經(jīng)用傳統(tǒng)的CAS方法來實(shí)現(xiàn)雙精度的原子操作，但是這種方法仍然會比CUDA運(yùn)行時(shí)環(huán)境提供的單精度原子方法慢得多。為了能擁有單精度運(yùn)算的快速和雙精度運(yùn)算的精確，本文將二者結(jié)合了起來。

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

Fig.5　Duplications in two dimensions圖5　在兩個(gè)維度上通過冗余解決數(shù)據(jù)沖突

5.2　粒子排序算法

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

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

Fig.6　Particle reordering圖6　粒子排序

5.3　CUDA-Aware MPI和GPUDirect RDMA

MPI能很好地和CUDA編程模型相兼容[13]。對于一般的MPI實(shí)現(xiàn)來說，只有指向host空間的指針可以被當(dāng)作參數(shù)傳遞給MPI函數(shù)。然而如果把MPI和CUDA相結(jié)合，就可以發(fā)送GPU的緩存，而不僅僅是CPU端的數(shù)據(jù)了。如果不利用CUDA-Aware的MPI實(shí)現(xiàn)，就需要利用cudaMemcpy函數(shù)把GPU內(nèi)存中的數(shù)據(jù)傳輸?shù)紺PU端的內(nèi)存。然而如果利用了這個(gè)技術(shù)，GPU端的緩存可以直接被MPI傳輸?shù)搅硪欢?。從下面的代碼對比中可以看出，使用了CUDA-Aware的MPI實(shí)現(xiàn)后，代碼會變得更簡潔，編程也更加容易。

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

（1）所有進(jìn)行數(shù)據(jù)傳輸?shù)牟僮鞫紩涣魉?/p>

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

遠(yuǎn)程DMA（RDMA）技術(shù)是在CUDA5.0中引入的GPUDirect技術(shù)中重要的一部分[15]，它連通了GPU和第三方應(yīng)用了PCI-E標(biāo)準(zhǔn)的硬件。GPU緩沖區(qū)中的數(shù)據(jù)可以不通過CPU端而直接被送到網(wǎng)卡進(jìn)行傳輸，從而消除了CPU到GPU和其他PCI-E設(shè)備的內(nèi)存帶寬占用。這樣也就顯著增大了GPU和其他節(jié)點(diǎn)的MPISendRecv的效率。在一個(gè)由許多GPU加速節(jié)點(diǎn)組成的集群中，嘗試應(yīng)用了這個(gè)技術(shù)來加速數(shù)據(jù)傳輸過程。通過實(shí)驗(yàn)發(fā)現(xiàn)，當(dāng)傳輸?shù)臄?shù)據(jù)量比較大時(shí)會有比較好的加速效果。隨著將來越來越多的加速卡和節(jié)點(diǎn)的投入使用，更好地利用這項(xiàng)技術(shù)來減少數(shù)據(jù)傳輸?shù)南膶兊迷絹碓街匾?/p>

6　實(shí)驗(yàn)結(jié)果和分析

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

Table 1　Software and hardware environment表1　實(shí)驗(yàn)的軟硬件環(huán)境

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

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

6.1　動態(tài)冗余算法

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

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

這里不考慮MPI通信時(shí)間隨著GPU數(shù)量增多而變化帶來的影響，因?yàn)槭褂萌哂鄡?yōu)化方法并不會改變通過MPI傳輸?shù)臄?shù)據(jù)數(shù)量。由于總的粒子數(shù)目是固定的，隨著使用GPU數(shù)的增加，每個(gè)GPU上的粒子數(shù)目會減少，這樣數(shù)據(jù)沖突的數(shù)量就會降低，冗余優(yōu)化的效果不突出，反而從冗余空間收集數(shù)據(jù)和更新數(shù)組的額外開銷變得明顯，導(dǎo)致加速比緩慢下降。

6.2　混合精度算法

總體上來說，隨著GPU數(shù)目的增多，混合精度計(jì)算方法帶來的加速收益會變得越來越明顯，當(dāng)GPU數(shù)量超過8后，加速效果會超過動態(tài)冗余算法。這是因?yàn)殡m然總的計(jì)算量是固定的，但是當(dāng)GPU數(shù)目較少時(shí)，每個(gè)GPU上的計(jì)算量就會較多，這樣進(jìn)行數(shù)據(jù)精度轉(zhuǎn)換帶來的額外開銷會比較大。當(dāng)每個(gè)GPU上的粒子數(shù)變少后，由于各個(gè)GPU上精度轉(zhuǎn)換函數(shù)是并行執(zhí)行的，精度轉(zhuǎn)換的總開銷也會相應(yīng)降低，加速效果就會變得越來越明顯。

6.3　RDMA

對于GPUDirect RDMA，本文通過改變MPI線程的數(shù)目來改變需要傳輸?shù)臄?shù)據(jù)量。發(fā)現(xiàn)使用GPUDirect RDMA可以顯著減少數(shù)據(jù)傳輸?shù)臅r(shí)間。如圖9所示，當(dāng)傳輸?shù)臄?shù)據(jù)量大于3 KB時(shí)，可以獲得2.8倍的加速比。但是當(dāng)傳輸?shù)臄?shù)據(jù)量為0.9 KB時(shí)，加速比只有1.37。這是因?yàn)楫?dāng)需要傳輸?shù)臄?shù)據(jù)量比較少時(shí)，RDMA的優(yōu)勢不明顯。待傳輸?shù)臄?shù)據(jù)量越大，越多的數(shù)據(jù)要傳輸?shù)搅硪粋€(gè)MPI域中。對于傳統(tǒng)的MPISendRecv方法來說，就需要更多的時(shí)間將數(shù)據(jù)從CPU端傳輸?shù)紾PU端。

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

這里需要指出的一點(diǎn)是，對于非RDMA版本，MPI通信時(shí)間的計(jì)算不僅僅包括純粹的數(shù)據(jù)傳輸時(shí)間，還包括收集數(shù)據(jù)，將數(shù)據(jù)放到數(shù)組中和把接收到的數(shù)據(jù)更新到正確的數(shù)組位置中的時(shí)間，因?yàn)镽DMA版本并不需要這些操作。

6.4　粒子排序

本文的粒子排序算法的中心思想是引入了一個(gè)參數(shù)empty_ratio，目的是在性能達(dá)到或者接近最優(yōu)的情況下找到這個(gè)參數(shù)的最小值。用8個(gè)MPI進(jìn)程來測試這個(gè)參數(shù)不同取值情況下的程序性能，結(jié)果如圖10所示。實(shí)驗(yàn)發(fā)現(xiàn)當(dāng)empty_ratio小于0.05時(shí)，性能受這個(gè)參數(shù)的影響很大。然而當(dāng)大于0.05時(shí)，性能接近于穩(wěn)定。因?yàn)檫@個(gè)參數(shù)越小代表著內(nèi)存空間的使用越少，所以可以通過設(shè)置這個(gè)參數(shù)為0.05來達(dá)到性能和內(nèi)存使用之間的平衡。

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

7　結(jié)論

本文研究了一種把整體的基于PIC方法的激光等離子體相互作用模擬的代碼移植到GPU端的可行方法，并且對比原始的CPU代碼取得了可觀的加速比?；谶@個(gè)初始的GPU版本，本文介紹了一系列的優(yōu)化方法來加速基于原子操作的包含大量數(shù)據(jù)沖突的Scatter階段，包括動態(tài)冗余算法、混合精度計(jì)算方法和一種參數(shù)化的粒子排序方法。也嘗試?yán)貌⑶以u估了GPUDirect RDMA方法在集群中加速M(fèi)PI通信時(shí)間的效果。發(fā)現(xiàn)當(dāng)數(shù)據(jù)傳輸量大于一定閾值的情況下，這種方法能顯著減少M(fèi)PI通信時(shí)間。相信這些優(yōu)化方法也能應(yīng)用于其他基于PIC方法的物理學(xué)模擬代碼，并且也對激光等離子體相互作用的研究發(fā)展有著非常重要的意義。

致謝林新華特別致謝日本學(xué)術(shù)振興會JSPS的RONPAKU項(xiàng)目資助。感謝NVIDIA GCOE的支持。感謝上海交通大學(xué)物理與天文系的翁蘇明教授的指導(dǎo)和幫助。

[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.