使用gpu并行實(shí)現(xiàn)的基于pic模型的加速器仿真方法
【專利摘要】本發(fā)明公開了使用GPU并行實(shí)現(xiàn)的基于PIC模型的加速器仿真方法，包括：將初始化信息從主機(jī)復(fù)制到計(jì)算節(jié)點(diǎn)的GPU中；根據(jù)初始化信息，確定粒子位置與網(wǎng)格的對應(yīng)關(guān)系；根據(jù)粒子位置與網(wǎng)格的對應(yīng)關(guān)系，計(jì)算每個網(wǎng)格內(nèi)所有粒子在網(wǎng)格上的電荷密度權(quán)重，得到網(wǎng)格的電荷密度分布；根據(jù)網(wǎng)格的電荷密度分布計(jì)算網(wǎng)格的電勢分布，并根據(jù)網(wǎng)格的電勢分布計(jì)算網(wǎng)格電場分布；計(jì)算每個粒子在電場作用下的運(yùn)動變化，并更新每個粒子的運(yùn)動狀態(tài)；以及用每個粒子的更新的運(yùn)動狀態(tài)代替初始化信息，迭代地執(zhí)行以上步驟，直到粒子運(yùn)動狀態(tài)滿足設(shè)計(jì)需求。根據(jù)本發(fā)明實(shí)施例，能夠解決現(xiàn)有PIC模型加速器仿真算法運(yùn)算速度低、成本高等技術(shù)問題。
【專利說明】使用GPU并行實(shí)現(xiàn)的基于PIC模型的加速器仿真方法
【技術(shù)領(lǐng)域】
[0001]本發(fā)明涉及加速器物理模擬仿真【技術(shù)領(lǐng)域】，具體涉及一種使用圖形處理單元(GPU)加速卡的計(jì)算機(jī)實(shí)現(xiàn)基于質(zhì)點(diǎn)網(wǎng)格(PIC)模型下對加速器中粒子加速過程模擬仿真方法。
【背景技術(shù)】
[0002]加速器模擬是加速器物理中一個非常重要的組成部分，細(xì)致而科學(xué)的加速器模擬對加速器設(shè)計(jì)和加速器調(diào)試起著至關(guān)重要的作用。由于在加速器中粒子數(shù)量巨大，粒子間相互作用緊密，以及加速器組成結(jié)構(gòu)復(fù)雜，都對加速器模擬技術(shù)提出了很高的要求。選取合適的模擬模型，采用高效的模擬算法，利用高性能的計(jì)算硬件，可以更加真實(shí)、細(xì)微和直觀地模擬出真實(shí)加速器環(huán)境中粒子的運(yùn)動形態(tài)。
[0003]目前，加速器模擬算法的實(shí)現(xiàn)以CPU為主，這些方法由于CPU計(jì)算能力不足導(dǎo)致計(jì)算規(guī)模不足，在模擬過程中，極大的限制了計(jì)算速度和模擬規(guī)模，在可以接受的機(jī)時(shí)內(nèi)只能計(jì)算很小的空間尺寸和時(shí)間尺寸。而基于NVIDIA公司開發(fā)的CUDA架構(gòu)的GPU并行編程為高性能運(yùn)算提供一種方便、直接又經(jīng)濟(jì)的解決方案。CUDA架構(gòu)的編程模式，依靠其高效的并行模式，方便的庫函數(shù)，為加速器模擬提供了新的契機(jī)。

【發(fā)明內(nèi)容】

[0004]本發(fā)明的目的在于，提供一種使用GPU加速卡并行實(shí)現(xiàn)基于PIC模型的加速器過程模擬的方法，以解決現(xiàn)有PIC模型加速器模擬算法運(yùn)算速度低、成本高等技術(shù)問題。
[0005]根據(jù)本發(fā)明一方面，提出了一種使用圖形處理單元(GPU)實(shí)現(xiàn)的基于質(zhì)點(diǎn)網(wǎng)格(Pic)模型的加速器仿真方法，包括:
[0006]a.在主機(jī)中產(chǎn)生初始化信息，并將初始化信息從主機(jī)復(fù)制到計(jì)算節(jié)點(diǎn)的GPU中，GPU包括多個流處理器；
[0007]在GPU中，利用多個流處理器并行執(zhí)行以下步驟:
[0008]b.根據(jù)初始化信息，確定粒子位置與網(wǎng)格的對應(yīng)關(guān)系；
[0009]c.根據(jù)粒子位置與網(wǎng)格的對應(yīng)關(guān)系，計(jì)算每個網(wǎng)格內(nèi)所有粒子在網(wǎng)格上的電荷密度權(quán)重，得到網(wǎng)格的電荷密度分布；
[0010]d.根據(jù)網(wǎng)格的電荷密度分布計(jì)算網(wǎng)格的電勢分布，并根據(jù)網(wǎng)格的電勢分布計(jì)算網(wǎng)格電場分布；
[0011]e.計(jì)算每個粒子在電場作用下的運(yùn)動變化，并更新每個粒子的運(yùn)動狀態(tài)；以及
[0012]f.用每個粒子的更新的運(yùn)動狀態(tài)代替初始化信息，迭代地執(zhí)行步驟b到e，直到粒子運(yùn)動狀態(tài)滿足設(shè)計(jì)需求。
[0013]在一個實(shí)施例中，在步驟b和e中按照一個流處理器對應(yīng)一個粒子的方式進(jìn)行GPU并行處理，并且在步驟c和d中按照一個流處理器對應(yīng)一個網(wǎng)格的方式進(jìn)行GPU并行處理。
[0014]在一個實(shí)施例中，初始化信息包括三維仿真空間劃分所得的網(wǎng)格數(shù)目、粒子的數(shù)目、粒子的三維位置和速度。
[0015]在一個實(shí)施例中，步驟b包括:確定每個粒子的位置所在的網(wǎng)格的編號并存儲在數(shù)組中；根據(jù)確定的編號對數(shù)組中的粒子位置進(jìn)行排序，使在同一個網(wǎng)格內(nèi)的所有粒子位置連續(xù)排列；在排序后的數(shù)組中獲取每個網(wǎng)絡(luò)內(nèi)粒子的開始位置和結(jié)束位置。
[0016]在一個實(shí)施例中，在步驟b中采用包括多個并行線程的線程塊，每個線程塊處理預(yù)定數(shù)目的網(wǎng)格，并且線程塊共享訪問GPU中的共享內(nèi)存。
[0017]在一個實(shí)施例中，在步驟d中，對電荷密度分布進(jìn)行三維傅立葉變換，根據(jù)頻域電荷密度分布計(jì)算網(wǎng)格的頻域電勢分布，并對頻域電勢分布進(jìn)行三維傅立葉逆變換，以得到網(wǎng)格的電勢分布。
[0018]在一個實(shí)施例中，步驟e包括:計(jì)算每個粒子在電場作用下的受力和加速度，并更新每個粒子的三維速度和位置。
[0019]在一個實(shí)施例中，步驟e還包括:更新所述數(shù)組中的粒子位置，對數(shù)組中所有粒子位置排序，并更新每個網(wǎng)格內(nèi)粒子的開始位置和結(jié)束位置。
[0020]在一個實(shí)施例中，利用排序后的數(shù)組，線程塊對連續(xù)排列的同一個網(wǎng)格內(nèi)的所有粒子位置進(jìn)行合并訪問。
[0021]在一個實(shí)施例中，在步驟d中利用紋理內(nèi)存綁定的方法，根據(jù)網(wǎng)格的電勢分布計(jì)算網(wǎng)格的電場分布。
[0022]在一個實(shí)施例中，在步驟e中，改變線程塊的大小，并利用改變后的線程塊的大小計(jì)算每個粒子的運(yùn)動變化。
[0023]本發(fā)明提供一種了使用GPU實(shí)現(xiàn)基于PIC模型的加速器模擬的方法，利用快速發(fā)展的GPU加速卡技術(shù)，采用例如NVIDIA公司的CUDA等并行計(jì)算架構(gòu)來并行實(shí)現(xiàn)該算法，可以更為真實(shí)地模擬加速器中束團(tuán)運(yùn)動過程，其實(shí)現(xiàn)效率高，成本低，可以更為便捷，快速地完成加速器模擬過程。此外，在算法結(jié)構(gòu)設(shè)計(jì)中結(jié)合GPU硬件，本發(fā)明設(shè)計(jì)出符合GPU并行模式的算法結(jié)構(gòu)，通過對粒子并行和網(wǎng)格并行的把握，充分利用GPU多線程和內(nèi)存共享的優(yōu)點(diǎn)，最大限度實(shí)現(xiàn)對PIC模型效率的提高。本發(fā)明有效地利用了現(xiàn)有基于CUDA架構(gòu)的CUFFT 庫(參見“CUDA FFT (CUFFT) Library DocumentationNVIDIA, 2012.)和 THRUST 庫(參見“CUDA TOOLKIT Documentation”，NVIDIA2012.)，采用粒子排序算法，在高效實(shí)現(xiàn)加速器模擬的同時(shí)，充分提高了程序的可移植性和可擴(kuò)展性。最后，本發(fā)明在實(shí)現(xiàn)過程中，利用共享內(nèi)存提高數(shù)據(jù)運(yùn)算速度、優(yōu)化數(shù)據(jù)結(jié)構(gòu)，利用符合GPU硬件架構(gòu)的合并訪問技術(shù)減少顯存讀取次數(shù)，合理利用紋理內(nèi)存，提高關(guān)鍵內(nèi)核函數(shù)的運(yùn)算速度，很大程度提高了整體程序的運(yùn)算效率。
【專利附圖】

【附圖說明】
[0024]下面通過附圖和實(shí)施例，對本發(fā)明的技術(shù)方案做進(jìn)一步的詳細(xì)描述。附圖中:
[0025]圖1是實(shí)施根據(jù)本發(fā)明實(shí)施例的使用GPU的基于PIC模型的加速器仿真方法的系統(tǒng)結(jié)構(gòu)的示意圖；以及
[0026]圖2是根據(jù)本發(fā)明實(shí)施例的使用GPU并行實(shí)現(xiàn)的基于PIC模型的加速器仿真方法的示意流程圖?！揪唧w實(shí)施方式】
[0027]以下結(jié)合附圖對本發(fā)明的優(yōu)選實(shí)施例進(jìn)行說明，應(yīng)當(dāng)理解，此處所描述的優(yōu)選實(shí)施例僅用于說明和解釋本發(fā)明，并不用于限定本發(fā)明。
[0028]圖1是可以用于實(shí)施根據(jù)本發(fā)明實(shí)施例的使用GPU的基于PIC模型的加速器仿真方法的系統(tǒng)結(jié)構(gòu)的示意圖。如圖1所示，該系統(tǒng)包括主機(jī)10、和計(jì)算節(jié)點(diǎn)20。雖然圖中未示出，但是可以理解，主機(jī)10可以包括用戶接口，例如鍵盤、觸摸板等，通過用戶接口，用戶可以輸入仿真所需的初始化信息、條件等。在圖1中，主機(jī)10和計(jì)算節(jié)點(diǎn)20之間經(jīng)由例如纜線等直接連接。然而，本發(fā)明實(shí)施例也可以采用其他任何適當(dāng)?shù)倪B接方式。例如，主機(jī)10可以經(jīng)由局域網(wǎng)或廣域網(wǎng)等與計(jì)算節(jié)點(diǎn)20通信，這使得能夠遠(yuǎn)程進(jìn)行仿真實(shí)驗(yàn)。例如，計(jì)算節(jié)點(diǎn)中的GPU可以通過PC1-E插口直接與主機(jī)的主板相連，完成CPU內(nèi)存和GPU顯存間的數(shù)據(jù)通信，并由該但GPU完成模擬過程。在仿真過程中，數(shù)據(jù)一旦從內(nèi)存中拷貝到GPU顯存中，所有計(jì)算均在GPU顯存中完成，CPU內(nèi)存可以只負(fù)責(zé)數(shù)據(jù)從顯存中提取回來和寫入磁盤當(dāng)中，或者進(jìn)行屏幕打印。可選地，圖1的系統(tǒng)還可以包括外部的存儲設(shè)備(未示出)，可以與主機(jī)連接，可以存儲例如各個計(jì)算節(jié)點(diǎn)的計(jì)算結(jié)果，以防止死機(jī)、斷電等意外情況發(fā)生導(dǎo)致數(shù)據(jù)丟失。
[0029]在一個實(shí)施例中，計(jì)算節(jié)點(diǎn)20可以是有GPU加速卡的高性能集群。在一個實(shí)施例中，計(jì)算節(jié)點(diǎn)20都具有GFllO核心以上的NVIDIA通用計(jì)算卡。計(jì)算節(jié)點(diǎn)20可以進(jìn)行基于NVIDIA公司開發(fā)的CUDA架構(gòu)的GPU并行編程。
[0030]在一個實(shí)施例中，主機(jī)10可以由例如包括中央處理單元(CPU)的通用計(jì)算機(jī)實(shí)現(xiàn)。
[0031]以上系統(tǒng)僅僅是本發(fā)明基本構(gòu)思的一種實(shí)現(xiàn)。本領(lǐng)域技術(shù)人員可以理解，上述各個部件的功能可以進(jìn)行再分配或組合，以形成其他的系統(tǒng)構(gòu)架。此外，如果功能足夠強(qiáng)大，上述各個部件的功能可以集成到單個計(jì)算機(jī)或工作站中。
[0032]下面結(jié)合圖1的系統(tǒng)結(jié)構(gòu)，參照圖2來描述根據(jù)本發(fā)明實(shí)施例的使用GPU實(shí)現(xiàn)的基于Pic模型的加速器仿真方法。圖2示出了該加速器仿真方法的示意流程圖。
[0033]如圖2所示，在步驟202，在主機(jī)10中產(chǎn)生初始化信息，并將初始化信息從主機(jī)復(fù)制到計(jì)算節(jié)點(diǎn)20的GPU中。具體地，初始化信息可以包括關(guān)于網(wǎng)格和粒子的信息，例如三維仿真空間劃分所得的網(wǎng)格數(shù)目、粒子的數(shù)目、粒子的三維位置和速度。在一個實(shí)施例中，初始化信息從主機(jī)內(nèi)存復(fù)制到計(jì)算節(jié)點(diǎn)中的GPU設(shè)備內(nèi)存中，粒子初始信息可以包括粒子數(shù)目ns、粒子三維方向位置信息pos_x, pos_y, pos_z和三維方向速度信息v_z, v_y, v_z。網(wǎng)格初始信息可以包括將仿真空間劃分為Nx*Ny*Nz個網(wǎng)格，Nx, Ny, Nz分別表示在x, y, z方向空間劃分的網(wǎng)格數(shù)目。
[0034]在步驟204，在GPU中，根據(jù)初始化信息，確定粒子位置與網(wǎng)格的對應(yīng)關(guān)系。計(jì)算節(jié)點(diǎn)的GPU對每個空間網(wǎng)格中的所有粒子信息進(jìn)行處理，其中計(jì)算節(jié)點(diǎn)GPU的每個流處理器負(fù)責(zé)處理一個對應(yīng)網(wǎng)格內(nèi)所有的粒子。需要讀取粒子的三維方向位置數(shù)組，即pos_X[]，pos_y[]，pos□，將粒子位置權(quán)重與空間網(wǎng)格格點(diǎn)進(jìn)行對應(yīng)。結(jié)合以上示例，步驟204具體包括如下過程:
[0035](I)首先在GPU中對粒子的三維方向位置口08_1,口08_7,口08_2做預(yù)處理,獲得粒子所在的網(wǎng)格信息。開辟大小為粒子數(shù)目ns的cell_count_list[]數(shù)組,cell_count_list []數(shù)組用來存放各個粒子所處在的空間網(wǎng)格的編號:
[0036]Pos—to—cell <<< ns/256, 256 > > > (cell—count—list，pos—x，pos_y, pos_z);
[0037]在一個實(shí)施例中，進(jìn)行cuda相關(guān)函數(shù)調(diào)用時(shí)，每一線程塊處理256個粒子，共有ns/256個線程塊。
[0038](2)對步驟(I)中獲得的粒子cell—count—list信息進(jìn)行連續(xù)排序，使在同一個網(wǎng)格內(nèi)的所有粒子在GPU設(shè)備端內(nèi)存空間連續(xù)排列。這一過程，釆用TRUST庫中的sort函數(shù)并行完成。THRUST庫函數(shù)實(shí)現(xiàn)排序可以有效利用GPU設(shè)備的流處理器并行計(jì)算，不需要對CPU與GPU端的數(shù)據(jù)進(jìn)行交換，避免了數(shù)據(jù)傳輸過程的時(shí)間開銷。
[0039]同時(shí)，在GPU中利用THRUST庫中的sort—by—key函數(shù)并行初始化排序索引cell—count—index[]，即將數(shù)組的第i個值賦值為i，i大于等于O小于等于ns。cell—count—index的意義是可以在排序后找回排序前粒子分布的原始信息:
[0040]Thrust::sort_by_key (cell_count_list，cell_count_list+ns，cell_count_index)。
[0041](3)在步驟⑵中得到cell—count—list排序的信息后，為了方便對同一網(wǎng)格內(nèi)的所有粒子統(tǒng)計(jì)，需要用search—start函數(shù)并行獲取每個網(wǎng)格內(nèi)粒子的開始位置cell—count—start []與結(jié)束位置cell—count—end[]，此時(shí)GPU使用ns/256個由256個線程組成的線程塊執(zhí)行這一過程，可以表示為:
[0042]Search—start <<< ns/256, 256 >>> (cell—count—list，cell—count—start，cell_count_end)；
[0043]此外，在上述獲取開始和結(jié)束位置開始前，使用cudaMemset函數(shù)重置cell—count—start []和 cell—count—end口全部為 0。
[0044]在Search—start 的 kernel 中，利用 GPU 中共享內(nèi)存(shared memory) shared—cell—count來提高計(jì)算效率，Shared—cell—count的大小定義為257個int型。
[0045]在GPU 硬件架構(gòu)中，共享內(nèi)存(參見 “NVIDIA CUDA C Programming Guide”，NVIDIA, 2012.)位于芯片上，所以共享內(nèi)存要比本地和全局內(nèi)存空間快得多。實(shí)際上，對于warp的所有線程，只要在線程之間沒有任何庫沖突，訪問共享內(nèi)存與訪問寄存器一樣快。所以，線程塊中可以共享訪問，訪問速度也遠(yuǎn)遠(yuǎn)高于全局顯存，可以極大地提高線程塊內(nèi)共享臨時(shí)數(shù)據(jù)的訪問速度，這里將每個線程中的cell—count—list對應(yīng)到shared—cell—count的共享顯存中，利用共享內(nèi)存的優(yōu)點(diǎn)，高效實(shí)現(xiàn)search—start操作。
[0046]Search—start的具體實(shí)現(xiàn)過程為:
[0047]每一個線程負(fù)責(zé)將第cell—count—list [blockldx.x*256+threadldx.x]的數(shù)值記錄到shared—cell—count [threadldx.x]上，而第0個線程還需同時(shí)負(fù)責(zé)將shared—cell—count [ (blockldx.x+1) *256]記錄到 shared—cell—count [256]的位置上。
[0048]此時(shí)，編號為blockldx.x*256+threadldx.x 的流處理器會將 shared—cell—count [threadldx.x]與 shared—cell—count [threadldx.x+1]兩個值進(jìn)行比較，若其值不一樣，則分別記錄為:
[0049]cel I_cοunt_start [blockldx.x*256 + threadldx.x] = shared_ce11_count[threadldx.x][0050]和
[0051]cell_count_start [blockldx.x*256+threadldx.x_l] = shared_cell_count[threadldx.x]
[0052](4)、在GPU中利用cell_count_index索引重新更新粒子的位置信息,仍以256個線程為一組線程塊，共ns/256個并行線程,編號為I = blockldx.x*256+threadldx.x的流處理器處理第i個粒子更新后的粒子三維位置和速度信息:
[0053]Pos_x_new[I] = pos_x[cell_count_index[i]]；
[0054]…
[0055]v_x_new[I] = pos_x[cell_count_index[i]]；
[0056]此時(shí)每個流處理器的內(nèi)核函數(shù)中包含cell_count_end[i] -cell_count_start [i]次循環(huán)。在如下描述的計(jì)算網(wǎng)格的電荷密度分布數(shù)組rho[]的過程中，利用更新過的粒子位置信息依次將編號從cell_count_start[i]至cell_count_end[i]的粒子的電荷密度權(quán)重存儲在電荷密度分布數(shù)組的相應(yīng)元素rho[i]中。
[0057]在CUDA架構(gòu)下，由GPU硬件架構(gòu)設(shè)計(jì)，全局顯存的讀入過程，是根據(jù)線程塊的大小批量讀入至內(nèi)核函數(shù)當(dāng)中的，數(shù)據(jù)在內(nèi).存中存儲時(shí)，如果線程對應(yīng)數(shù)據(jù)在內(nèi)存中連續(xù)時(shí)，可以減少內(nèi)存數(shù)據(jù)的載入次數(shù)，這個過程叫做合并訪問。在這種算法設(shè)計(jì)下的，粒子在排序后，存儲順序由粒子的網(wǎng)格編號決定，同一個網(wǎng)格編號的粒子內(nèi)存連續(xù)，在內(nèi)存讀取過程中，達(dá)成內(nèi)存合并訪問，可以有效的減少線程塊block中的內(nèi)存讀取次數(shù)，極大的提高了內(nèi)存訪問速度，有效的提高了 GPU并行權(quán)重的效率。
[0058]在步驟206，根據(jù)粒子位置與網(wǎng)格的對應(yīng)關(guān)系，計(jì)算每個網(wǎng)格內(nèi)所有粒子在網(wǎng)格上的電荷密度權(quán)重，得到網(wǎng)格的電荷密度分布。具體地，初始化電荷密度數(shù)組rho[]，利用cudaMemset函數(shù)將初值全部賦值為O。GPU的每個線程塊包含256個線程，共分配(Nx*Ny*Nz/256)個線程塊，每個流處理器I = blockldx.x*256+threadldx.x處理編號為i的網(wǎng)格內(nèi)所有粒子的電荷密度權(quán)重。
[0059]在步驟208，根據(jù)網(wǎng)格的電荷密度分布計(jì)算網(wǎng)格的電勢分布，并根據(jù)網(wǎng)格的電勢分布計(jì)算網(wǎng)格電場分布。在GPU中，根據(jù)空間網(wǎng)格的電荷密度分布，使用CUFFT庫函數(shù)，利用GI3U流處理器對應(yīng)空間網(wǎng)格的并行方式并行地將電荷密度分布轉(zhuǎn)換為頻域電荷密度分布，通過頻域電荷密度分布求解空間網(wǎng)格格點(diǎn)處的頻域電勢分布，再使用CUFFT庫函數(shù)，將頻域電勢分布轉(zhuǎn)換為網(wǎng)格點(diǎn)上的電勢分布，進(jìn)而并行求解出空間網(wǎng)格的電場分布。步驟208具體包括如下過程。
[0060](I)在GPU設(shè)備端，利用cufft庫函數(shù)，對電荷密度分布作三維傅立葉正變換，變換后，得到頻域電荷密度分布。需注意，在變換之前，應(yīng)將電荷密度分布格式由float*轉(zhuǎn)換為cuffiComplex*,,使之成為可以使用CUFFT庫函數(shù)的數(shù)據(jù)格式。
[0061]可以采用如下的操作來實(shí)現(xiàn)傅立葉變化，
[0062]cufftPlan3d(&plan, Nz, Ny, Ny, CUFFT_C2C)；
[0063]cufftExecC2C (plan, rho, rho_fft,，CUFFT_F0RWARD)；
[0064]在這個變換過程中，Nz, Ny, Nx的順序不得改變，它是由三維網(wǎng)格向一維網(wǎng)格轉(zhuǎn)換過程的順序決定的；變換過程中，輸入數(shù)據(jù)為rho，輸出數(shù)據(jù)為rho_fft ;變換標(biāo)識為CUFFT_FORWARD,也就是采用傅立葉正變換。[0065]變換后，得到頻域的電荷密度分布rho_fft []；
[0066](2)根據(jù)步驟⑴中獲得的頻域電荷密度，可在GPU設(shè)備中并行地求解頻域電勢分布，即:
[0067]Rho_to_phi <<< (Nx*Ny*Nz/256), 256 >>> (phi_fit, rho_fft)；
[0068]其中，phi_fft為頻域電勢分布，rho_fft為頻域電荷密度。在內(nèi)核函數(shù)中，流處理I = blockldx.x*256+threadldx.x被用于求解第i個網(wǎng)格內(nèi)頻域電勢，通過計(jì)算即可得到空間內(nèi)所有網(wǎng)格格點(diǎn)上的頻域電勢分布。
[0069]該過程需要大量使用GPU中的臨時(shí)寄存器，為了最高效地發(fā)揮GPU的處理能力，采用每個線程塊中包含256個線程的線程塊劃分方式。
[0070](3)根據(jù)步驟(2)得到的頻域電勢分布，繼續(xù)利用CUFFT庫函數(shù)，對頻域電勢分布做傅立葉逆變換:
[0071]cufftPlan3d(&plan, Nz, Ny, Ny, CUFFT_C2C)；
[0072]cufftExecC2C(plan, phi_fit, phi, CUFFT_INVERSE)；
[0073]在這個變換過程中，Nz, Ny, Nx的順序不得改變，它是由三維網(wǎng)格向一維網(wǎng)格轉(zhuǎn)換過程的順序決定的；變換過程中，輸入數(shù)據(jù)為phi_fft，輸出數(shù)據(jù)為phi ;變換標(biāo)識為CUFFT_INVERSE,也就是采用傅立葉逆變換。
[0074]變換后，需要對網(wǎng)格點(diǎn)上的電勢分布進(jìn)行數(shù)據(jù)格式的轉(zhuǎn)化，即將cufftComplex*轉(zhuǎn)化為float*，轉(zhuǎn)化后，得到電勢分布。
[0075](4)分配(Nx*Ny*Nz/256)個線程塊，每塊線程塊由256個線程組成，將步驟(3)中得到的空間電勢分布分別在X，y和z方向上求解空間網(wǎng)格電場分布。
[0076]在這個過程中，采用紋理內(nèi)存(參見“NVIDIA CUDA C Programming Guide”,NVIDIA, 2012)綁定的方法實(shí)現(xiàn)。紋理內(nèi)存空間具有高速緩存，所以紋理拾取僅在高速緩存缺失時(shí)，耗費(fèi)從設(shè)備內(nèi)存中的一個內(nèi)存讀取，否則它僅耗費(fèi)從紋理高速緩存中的一個讀取。紋理高速緩存針對2D空間局部性進(jìn)行了優(yōu)化，所以讀取緊密相鄰的紋理地址的同一 WARP的線程將達(dá)到最佳性能。此外，它還設(shè)計(jì)用于流水化具有恒定延遲的拾取。紋理內(nèi)存綁定，有更高的訪問帶寬，而且不受訪問模式的約束，在更快尋址的同時(shí)，可以隱藏計(jì)算時(shí)間，有時(shí)候會改善應(yīng)用程序之行隨機(jī)訪問數(shù)據(jù)的性能。結(jié)合上述示例，紋理內(nèi)存綁定的方法可以包括如下步驟:
[0077]a)綁定電勢分布為紋理內(nèi)存:
[0078]cudaBindTexture(O, rt, phi, Nx*Ny*Nz)；
[0079]b)利用電勢分布phi □,求出空間中的電場分布Ex[], Ey [], Ez []:
[0080]phi_to_Ex <<< Nx*Ny*Nz/256, 256 >>> (phi, Ex)；
[0081]...[0082]phi_to_Ez <<< Nx*Ny*Nz/256, 256 >>> (phi, Ez)；
[0083]在計(jì)算過程中，一共有Nx*Ny*Nz個格點(diǎn),流處理對應(yīng)網(wǎng)格,分別在x、y和z方向計(jì)算每個格點(diǎn)處的電場分布，GPU內(nèi)核劃分為Nx*Ny*Nz/256個線程塊，每個線程塊對應(yīng)256個線程。在三維方向計(jì)算電場分布的時(shí)候，由于電勢Phi []是以一維數(shù)組方式存儲，電勢phi □在y方向和z方向分布不連續(xù),紋理內(nèi)存的使用對Y方向和z方向的電場分布Ey和Ez的計(jì)算會有較大提高。[0084]c)解除紋理內(nèi)存綁定:
[0085]cudaUnBindTexture (rt)；
[0086]在使用過紋理內(nèi)存后，解除phi 口的紋理內(nèi)存綁定。
[0087]在步驟210，計(jì)算每個粒子在電場作用下的運(yùn)動變化，并更新每個粒子的運(yùn)動狀態(tài)。具體地，根據(jù)步驟208中所求的電場分布，如下并行計(jì)算每個粒子在所受電場力作用下的位置和速度變化，并更新粒子位置及速度變化:
[0088]Vel_ < < < ns/128,128 > > > (v_x, v_y, v_z, Ex, Ey, Ez, pos_x, pos_y, pos_z,v_x_new, v_y_new, v_z_new)；
[0089]Pos <<< ns/256, 256 >>> (、pos_x, pos_y, pos_z, v_x, v_y, v_z, pos_x_new, pos_y_new, pos_z_new)；
[0090]其中 pos_x_new, pos_y_new, pos_z_new, v_x_new, v_y_new, v_z_new 為步驟 204中更新后得到的粒子位置及粒子速度，pos_x, pos_y, pos_z, v_x, v_y, v_z是此次計(jì)算后更新得到的粒子位置及速度。
[0091 ] 在CUDA架構(gòu)中，由于GPU硬件的寄存器數(shù)量有限，在內(nèi)核函數(shù)中，如果臨時(shí)變量過多，則需要減少線程塊中的線程數(shù)目。若同一個線程塊中的線程數(shù)目過大，則會造成GPU中臨時(shí)寄存器數(shù)目不夠，使得計(jì)算速度急速下降甚至計(jì)算出錯。
[0092]在計(jì)算速度變化的內(nèi)核函數(shù)中，需要利用GPU中較多的寄存器存儲臨時(shí)變量，需要控制線程塊中的線程數(shù)目，所以每個線程塊大小劃分為128個線程，線程塊個數(shù)為ns/128 個。
[0093]在步驟212，判斷仿真結(jié)果是否滿足設(shè)計(jì)需求。如果不滿足，則將步驟210中更新后的粒子位置及速度變化值帶回到步驟204中繼續(xù)計(jì)算，直到仿真結(jié)果滿足設(shè)計(jì)需求。如果滿足，則以上過程結(jié)束。
[0094]可以將GPU設(shè)備內(nèi)存中的結(jié)果數(shù)據(jù)拷貝回主機(jī)內(nèi)存，并釋放相應(yīng)的CPU和GPU端內(nèi)存。
[0095]本發(fā)明實(shí)施例根據(jù)GPU多輕量計(jì)算核心的特點(diǎn)，對算法結(jié)構(gòu)及模型結(jié)構(gòu)進(jìn)行合理的安排與設(shè)計(jì)，將GPU中流處理器與算法中的網(wǎng)格、粒子有效地結(jié)合在一起，使PIC模型更加適合GPU并行模式，同時(shí)，充分提高了 GPU線程的使用效率，極大地減少了模擬過程的計(jì)算時(shí)間。
[0096]此外，本發(fā)明中，通過合理利用高速訪問的GPU中的共享內(nèi)存，優(yōu)化數(shù)據(jù)結(jié)構(gòu)，采用符合GPU硬件結(jié)構(gòu)的合并訪問，有效的減少顯存讀取次數(shù)，合理利用紋理內(nèi)存，對關(guān)鍵內(nèi)核函數(shù)提高30%左右的運(yùn)算效率，很大限度的發(fā)揮GPU計(jì)算的優(yōu)點(diǎn)。
[0097]本領(lǐng)域的技術(shù)人員可以對本發(fā)明進(jìn)行各種改動和變型而不脫離本發(fā)明的精神和范圍。這樣，倘若本發(fā)明的這些修改和變型屬于本發(fā)明權(quán)利要求及其等同技術(shù)的范圍之內(nèi)，則本發(fā)明也意圖包含這些改動和變型在內(nèi)。
【權(quán)利要求】
1.一種使用圖形處理單元GPU實(shí)現(xiàn)的基于質(zhì)點(diǎn)網(wǎng)格PIC模型的加速器仿真方法，包括: a.在主機(jī)中產(chǎn)生初始化信息，并將初始化信息從主機(jī)復(fù)制到計(jì)算節(jié)點(diǎn)的GPU中，GPU包括多個流處理器； 在GPU中利用多個流處理器并行執(zhí)行以下步驟: b.根據(jù)初始化信息，確定粒子位置與網(wǎng)格的對應(yīng)關(guān)系； c.根據(jù)粒子位置與網(wǎng)格的對應(yīng)關(guān)系，計(jì)算每個網(wǎng)格內(nèi)所有粒子在網(wǎng)格上的電荷密度權(quán)重，得到網(wǎng)格的電荷密度分布； d.根據(jù)網(wǎng)格的電荷密度分布計(jì)算網(wǎng)格的電勢分布，并根據(jù)網(wǎng)格的電勢分布計(jì)算網(wǎng)格電場分布； e.計(jì)算每個粒子在電場作用下的運(yùn)動變化，并更新每個粒子的運(yùn)動狀態(tài)；以及 f.用每個粒子的更新的運(yùn)動狀態(tài)代替初始化信息，迭代地執(zhí)行步驟b到e，直到粒子運(yùn)動狀態(tài)滿足設(shè)計(jì)需求。
2.根據(jù)權(quán)利要求1所述的方法，其中在步驟b和e中按照一個流處理器對應(yīng)一個粒子的方式進(jìn)行GPU并行處理，并且在步驟c和d中按照一個流處理器對應(yīng)一個網(wǎng)格的方式進(jìn)行GPU并行處理。
3.根據(jù)權(quán)利要求1或2所述的方法，其中初始化信息包括三維仿真空間劃分所得的網(wǎng)格數(shù)目、粒子的數(shù)目、粒子的三維位置和速度。
4.根據(jù)權(quán)利要求3所述的方法，其中步驟b包括: 確定每個粒子的位置所在的網(wǎng)格的編號并存儲在數(shù)組中； 根據(jù)確定的編號對數(shù)組中的粒子位置進(jìn)行排序，使在同一個網(wǎng)格內(nèi)的所有粒子位置連續(xù)排列； 在排序后的數(shù)組中獲取每個網(wǎng)絡(luò)內(nèi)粒子的開始位置和結(jié)束位置。
5.根據(jù)權(quán)利要求4所述的方法，其中在步驟b中采用包括多個并行線程的線程塊，每個線程塊處理預(yù)定數(shù)目的網(wǎng)格，并且線程塊共享訪問GPU中的共享內(nèi)存。
6.根據(jù)權(quán)利要求1或2所述的方法，其中，在步驟d中，對電荷密度分布進(jìn)行三維傅立葉變換，根據(jù)頻域電荷密度分布計(jì)算網(wǎng)格的頻域電勢分布，并對頻域電勢分布進(jìn)行三維傅立葉逆變換，以得到網(wǎng)格的電勢分布。
7.根據(jù)權(quán)利要求1或2所述的方法，其中步驟e包括:計(jì)算每個粒子在電場作用下的受力和加速度，并更新每個粒子的三維速度和位置。
8.根據(jù)權(quán)利要求4中所述的方法，其中步驟e還包括:更新所述數(shù)組中的粒子位置，對數(shù)組中所有粒子位置排序，并更新每個網(wǎng)格內(nèi)粒子的開始位置和結(jié)束位置。
9.根據(jù)權(quán)利要求4所述的方法，其中利用排序后的數(shù)組，線程塊對連續(xù)排列的同一個網(wǎng)格內(nèi)的所有粒子位置進(jìn)行合并訪問。
10.根據(jù)權(quán)利要求1或2所述的方法，其中在步驟d中利用紋理內(nèi)存綁定的方法，根據(jù)網(wǎng)格的電勢分布計(jì)算網(wǎng)格的電場分布。
11.根據(jù)權(quán)利要求1或2所述的方法，其中在步驟e中，改變線程塊的大小，并利用改變后的線程塊的大小計(jì)算每個粒子的運(yùn)動變化。
【文檔編號】G06F9/455GK103440163SQ201310413539
【公開日】2013年12月11日 申請日期:2013年9月9日 優(yōu)先權(quán)日:2013年9月9日
【發(fā)明者】楊磊, 張智磊, 李超, 齊新, 高笑菲 申請人:中國科學(xué)院近代物理研究所