魏曉輝, 朱 彤, 李洪亮, 李維山, 許天福

(1. 吉林大學(xué) 計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院, 長(zhǎng)春 130012; 2. 吉林大學(xué) 環(huán)境與資源學(xué)院, 長(zhǎng)春 130012)

隨著高性能并行計(jì)算技術(shù)的快速發(fā)展, 將地下多相流動(dòng)數(shù)值模擬與并行計(jì)算相結(jié)合, 提高地下多相流動(dòng)數(shù)值模擬能力, 已成為解決各種水文地質(zhì)與環(huán)境地質(zhì)問題的有效工具. 由美國勞倫斯伯克利實(shí)驗(yàn)室開發(fā)的TOUGHREACT[1]是當(dāng)前廣泛使用的解決地下多相流體運(yùn)動(dòng)與地球化學(xué)反應(yīng)運(yùn)移耦合問題的模擬程序. 但該程序在解決涉及大尺度、 高精度和高復(fù)雜性的大規(guī)模數(shù)值模擬問題時(shí), 應(yīng)用效果欠佳. 因此, 通過高性能并行計(jì)算技術(shù)加速TOUGHREACT的數(shù)值模擬過程意義重大, 其中一種方法就是利用圖形處理器(GPU)通用計(jì)算技術(shù)實(shí)現(xiàn)[2].

目前, 有關(guān)地下水流數(shù)值模擬的并行化研究主要基于集群系統(tǒng), 結(jié)合MPI技術(shù)和現(xiàn)有的并行線性求解軟件包(AZTEC, Petsc等)實(shí)現(xiàn). 如: Zhang等[3]利用TOUGH2-MP模擬了位于NevadaYucca山脈非飽和區(qū)域的水氣運(yùn)輸及熱傳導(dǎo); Huang等[4]結(jié)合MPI消息傳遞編程模型及Schwarz類區(qū)域分解法實(shí)現(xiàn)了MODFLOW與RT3D的并行化; Dong等[5]針對(duì)MOLDFLOW中的PCG(preconditioned conjugate gradient)軟件包, 基于多核CPU平臺(tái), 應(yīng)用OpenMP并行編程方法, 對(duì)修正的不完全Cholesky共軛梯度法(MICCG )進(jìn)行了部分并行化, 對(duì)多項(xiàng)式共軛梯度法(POLCG)進(jìn)行全部并行化, 實(shí)現(xiàn)了MODFLOW在共享存儲(chǔ)環(huán)境下的并行計(jì)算, 在8核心的并行計(jì)算機(jī)可以獲得1.40～5.31倍的加速比; Ji等[6]實(shí)現(xiàn)了一種基于CUDA的大規(guī)模地下水流模擬求解器, 對(duì)模擬過程中使用Jacobi預(yù)條件子的共軛梯度法進(jìn)行并行化, 并在不同的GPU上執(zhí)行, 相對(duì)于CPU串行程序達(dá)到了1.8～3.7倍的加速比.

本文通過對(duì)TOUGHREACT軟件的相關(guān)原理與并行性分析, 提出并實(shí)現(xiàn)了一種基于CUDA平臺(tái)的線性方程組并行求解策略, 在保證算法正確性和求解精度的前提下, 相比于串行程序取得了良好的加速效果. 然后將并行程序整合入TOUGHREACT軟件包中替代原有的求解器, 實(shí)現(xiàn)了基于CPU+GPU異構(gòu)體系的地下多相流動(dòng)數(shù)值模擬.

1 預(yù)備知識(shí)

1.1 GPU通用計(jì)算

基于GPU通用計(jì)算技術(shù)(general-purpose computing on graphics processing units, GPGPU)的并行策略與基于OpenMP的方式類似, 也是利用多線程實(shí)現(xiàn)對(duì)求解過程中費(fèi)時(shí)操作的并行計(jì)算, 而且GPU能開辟許多線程進(jìn)行計(jì)算, 更適合完成具有高度并行化的計(jì)算密集型任務(wù). 在GPGPU體系結(jié)構(gòu)下, 執(zhí)行并行程序時(shí)會(huì)產(chǎn)生CPU與GPU之間的數(shù)據(jù)通信. 因此, 進(jìn)行基于GPU的并行程序設(shè)計(jì)時(shí), 在最大限度挖掘程序并發(fā)性的同時(shí), 還要考慮如何控制CPU與GPU之間的通信開銷.

CUDA(compute unified device architecture)編程平臺(tái)以C語言為基礎(chǔ), 編程人員可使用熟悉的C語言調(diào)用CUDA中的庫函數(shù)完成GPU編程. 在CUDA程序架構(gòu)中, 程序執(zhí)行區(qū)域分為兩部分： Host與Device. Host指CPU, 在CPU上執(zhí)行的代碼稱為Host代碼； Device指GPU, 在GPU上執(zhí)行的代碼稱為Kernel代碼. CPU控制GPU執(zhí)行, 調(diào)度分配任務(wù), 而大量需要并行計(jì)算的工作則由GPU實(shí)現(xiàn). 圖1為CUDA程序的Kernel執(zhí)行與線程組織模型. 由圖1可見, 每個(gè)block內(nèi)的所有線程必須存在于同一個(gè)處理器核心中, 且共享該核心有限的存儲(chǔ)器資源, 因此, 一個(gè)block內(nèi)的線程數(shù)是有限的. 在目前的GPU上, 一個(gè)block可以包含多達(dá)1 024個(gè)線程. 實(shí)際執(zhí)行時(shí), 每個(gè)block內(nèi)的線程以32為單位組織成線程束(warp), warp是GPU上處理器資源調(diào)度的基本單元, 同一warp內(nèi)的線程天然同步.

圖1 CUDA程序Kernel執(zhí)行與線程組織模型Fig.1 CUDA program kernels and threads structure

1.2 空間離散及相關(guān)求解方法

TOUGHREACT采用積分有限差分法(integral finite difference method, IFDM)[7]將模擬區(qū)域離散成任意形狀的多面體. 計(jì)算過程中, 只需要單元的體積與面積及單元中心到各個(gè)面的垂直距離, 使得該方法在處理任意形狀的單元時(shí)不必考慮總體坐標(biāo)系統(tǒng), 同時(shí)也不受單元塊鄰近單元數(shù)的限制. 基于多相流具有很強(qiáng)的非線性特征, 程序利用Newton-Raphson迭代法求解非線性方程組, 并在計(jì)算過程中, 根據(jù)收斂速度調(diào)整時(shí)間步長(zhǎng). Newton迭代法的求解過程是一種數(shù)值逼近的計(jì)算方法, 迭代過程中, 每步迭代都要求解一個(gè)大型稀疏線性方程組. 共軛梯度類算法是目前的主流算法, 該方法需要的控制參量少, 迭代算法穩(wěn)定, 不易使解發(fā)散, 對(duì)不同類型的問題有較強(qiáng)的適應(yīng)性, 結(jié)合預(yù)處理技術(shù), 容易形成穩(wěn)定、 高效的線性方程求解器. TOUGHREACT提供給用戶3種預(yù)處理共扼梯度法求解稀疏矩陣線性方程組, 分別是雙共扼梯度求解、 Lanczos類的雙共扼梯度求解和穩(wěn)定雙共扼梯度求解.

2 數(shù)值模擬過程并行性分析

2.1 模擬過程分析及并行化

TOUGHREACT經(jīng)過不斷的完善擴(kuò)充, 源代碼十分龐雜, 移植工作量巨大. 因此, 本文采取將一部分熱點(diǎn)代碼優(yōu)先向GPU移植的策略. 先對(duì)程序結(jié)構(gòu)和執(zhí)行過程做簡(jiǎn)要介紹. 程序的數(shù)值模擬過程大致可分為兩部分： 每一個(gè)時(shí)間步, 先進(jìn)行地下多相流體運(yùn)動(dòng)的模擬計(jì)算并更新相應(yīng)網(wǎng)格的物理參數(shù), 然后進(jìn)行地球化學(xué)運(yùn)移計(jì)算, 二者的計(jì)算過程都涉及到數(shù)值方程求解, 但方程的組織形式和規(guī)模不同. 地球化學(xué)運(yùn)移過程中所形成的方程組是對(duì)每個(gè)網(wǎng)格計(jì)算單元進(jìn)行組建, 而地下多相流的方程組是在所有網(wǎng)格計(jì)算單元的層次上組建, 在數(shù)值計(jì)算方法上有差別, 本文主要對(duì)地下多相流體運(yùn)動(dòng)的數(shù)值模擬計(jì)算問題在GPU上并行實(shí)現(xiàn).

圖2 TOUGH2模塊化結(jié)構(gòu)Fig.2 TOUGH2 program module structure

TOUGHREACT沿用了TOUGH2[8]中求解多相流問題的相關(guān)模塊. TOUGH2的程序結(jié)構(gòu)設(shè)置為模塊化如圖2所示, 主要的水流和運(yùn)移模塊可以與不同水流屬性模型相互作用, 處理不同的多組分、 多相流系統(tǒng).

根據(jù)文獻(xiàn)[9]的分析, 程序的主要計(jì)算時(shí)間集中在對(duì)非線性方程組的Newton迭代過程(80%以上). 在此過程中, 組建線性方程組Ax=b的Jacobian系數(shù)矩陣A的時(shí)間約占30%, 線性方程組求解時(shí)間約占60%. 因此, 本文對(duì)其最耗時(shí)部分(稀疏線性方程組的求解)進(jìn)行GPU并行優(yōu)化. 圖3為程序執(zhí)行流程和基本任務(wù)劃分模塊.

圖3 主要模塊GPU計(jì)算流程Fig.3 Implementation of computing procedure on GPU

2.2 預(yù)處理穩(wěn)定雙共軛梯度法并行性分析

為了達(dá)到更好的收斂效果, 本文使用預(yù)處理穩(wěn)定雙共軛梯度法求解線性方程組時(shí), 先通過不完全LU分解(ILU)[10]對(duì)方程組的系數(shù)矩陣進(jìn)行預(yù)處理. 目前對(duì)稀疏矩陣ILU分解的并行化瓶頸主要在于現(xiàn)有算法的數(shù)據(jù)依賴性很強(qiáng), 導(dǎo)致可利用的并行性不足, 不適合在GPU上實(shí)現(xiàn). 本文中的不完全LU分解操作由CPU串行執(zhí)行.

下面對(duì)程序中涉及到求解線性方程組的穩(wěn)定雙共軛梯度法(Bi-conjugate gradient stabilized, BiCGSTAB)[10]進(jìn)行分析. 算法的描述參見文獻(xiàn)[10], 算法從初始解x(0)開始每步迭代都要進(jìn)行如下操作：

設(shè)系數(shù)矩陣A為n階可逆矩陣, 平均每行非零元個(gè)數(shù)為k,M(M=L×D×U)為通過不完全LU分解得到的預(yù)條件子, 其中:L和U分別是m階上、 下三角矩陣, 平均每行非零元數(shù)為s;D為m階對(duì)角矩陣. 每次迭代中各種計(jì)算所需次數(shù)及復(fù)雜度列于表1.

表1 穩(wěn)定的雙共軛梯度法中各計(jì)算步的計(jì)算復(fù)雜度分析

由表1可見, 綜合考慮各計(jì)算步驟的計(jì)算次數(shù)和復(fù)雜度, 稀疏矩陣與向量的乘積運(yùn)算(sparsematrix-vectormultiplication,SPMV), 向量?jī)?nèi)積和三角方程組求解是算法中最耗時(shí)部分, 對(duì)使用穩(wěn)定雙共軛梯度法求解的計(jì)算性能有很大影響.SPMV操作具有行無關(guān)的特性, 因此非常適合在GPU上并行實(shí)現(xiàn). 考慮到無論求解上三角方程組, 還是下三角方程組, 在回代過程中, 并行度會(huì)越來越小, 從m-1到1逐次遞減. 同時(shí)對(duì)于稀疏矩陣, 非零元素所占比例非常小, 并行度會(huì)更小. 因此, 本文在CPU上進(jìn)行三角方程組的求解, 稀疏矩陣與向量乘積(SPMV)等計(jì)算步驟均由GPU執(zhí)行. 此外, 由于向量?jī)?nèi)積運(yùn)算在GPU中運(yùn)行比SPMV遇到的性能瓶頸少很多, 所以優(yōu)化工作將集中在SPMV上.

3 穩(wěn)定雙共軛梯度法的GPU并行化

由上述分析可對(duì)并行程序中CPU和GPU的計(jì)算任務(wù)做出劃分, 將整個(gè)求解過程分成稀疏矩陣向量乘積運(yùn)算、 向量?jī)?nèi)積、 向量加減、 向量乘標(biāo)量和預(yù)處理方程組求解等幾個(gè)步驟. 為了獲得較好的整體優(yōu)化效果, 本文對(duì)不同操作采用不同策略.

3.1 稀疏矩陣存儲(chǔ)格式

對(duì)于線性代數(shù)方程組Ax=b, 其系數(shù)矩陣A是一個(gè)對(duì)稱正定的大型稀疏矩陣. 選擇合理的稀疏矩陣存儲(chǔ)格式對(duì)迭代求解計(jì)算性能的提高有顯著作用. 對(duì)稀疏矩陣存儲(chǔ)格式的定義應(yīng)包括非零元素的值和非零元素在矩陣中的位置信息. 常用的稀疏矩陣存儲(chǔ)格式是行壓縮(compressed sparse row, CSR)格式, CSR格式把一個(gè)二維矩陣壓縮到一維數(shù)組中, 對(duì)某個(gè)稀疏矩陣, 可以把全部非零元素記錄在數(shù)組Data中, 然后用數(shù)組Ptr記錄每行上第一個(gè)非零元素在數(shù)組Data中的索引, 每個(gè)非零元素對(duì)應(yīng)的列坐標(biāo)記錄在Indices數(shù)組中. 基于CSR格式的SPMV具有高度的并行性, 有利于使用GPU對(duì)SPMV進(jìn)行并行加速. 稀疏矩陣CSR存儲(chǔ)格式為

3.2 并行歸約求和

并行歸約(reduction)是一種可以實(shí)現(xiàn)并行求和、 求乘、 求最大(小)值等操作的方法, 用途廣泛. 本文使用并行歸約的方法對(duì)同一block(warp)內(nèi)線程求得的中間結(jié)果進(jìn)行并行求和. 初始階段每個(gè)線程讀入一個(gè)數(shù)據(jù), 在后續(xù)循環(huán)中, 每輪使用上一輪循環(huán)中一半的線程進(jìn)行tid與tid+s的求和, 其中s為線程跨度. 以half-warp(16個(gè)線程)為例, 有：

1)s=8. 由tid=0,1,2,3,4,5,6,7的線程執(zhí)行data[tid]+data[tid+8]計(jì)算;

2)s=4. 由tid=0,1,2,3的線程執(zhí)行data[tid]+data[tid+4]計(jì)算;

3)s=2. 由tid=0,1的線程執(zhí)行data[tid]+data[tid+2]計(jì)算;

4)s=1. 由tid=0的線程執(zhí)行data[tid]+data[tid+1]計(jì)算, 此時(shí)data[0]中保存的即為該half-warp輸入的16個(gè)數(shù)據(jù)之和.

圖4 規(guī)約求和Fig.4 Reduction procedure

由于下一輪循環(huán)中的線程要用到上一輪中的計(jì)算結(jié)果, 所以在每次循環(huán)中都要進(jìn)行一次線程同步操作, reduction的執(zhí)行過程如圖4所示.

3.3 稀疏矩陣向量乘并行化

本文借鑒文獻(xiàn)[11]的CSR格式稀疏矩陣與向量乘積在GPU中的實(shí)現(xiàn)方法, 利用warp內(nèi)線程天然同步特性, 由每個(gè)warp計(jì)算結(jié)果向量的一個(gè)元素, warp中的每個(gè)線程負(fù)責(zé)計(jì)算每個(gè)非零元與相應(yīng)標(biāo)量的乘積, 計(jì)算過程相互獨(dú)立, 可較大限度利用GPU處理器資源.

基于CSR格式的SPMV的Kernel.

thread_id =blockDim.x×blockIdx.x+threadIdx.x;； //每個(gè)線程的編號(hào)

warp_id=thread_id/32; //每個(gè)線程所在warp的編號(hào)

lane=thread_id&(32-1); //每個(gè)線程在當(dāng)前warp中的編號(hào)

_shared_vals[128]; //定義共享內(nèi)存數(shù)組;

row=warp_id; //每個(gè)warp計(jì)算結(jié)果向量的一個(gè)結(jié)果;

If (row

vals [threadIdx.x]=0;

for (jj=ptr [row]+lane;jj

vals[threadIdx.x]+=data[jj]×x[indices [jj]]; //warp內(nèi)線程同時(shí)計(jì)算非零元素與相 應(yīng)標(biāo)量的乘積, 中間結(jié)果存入數(shù)組vals中；

for(0≤lane<32);

Reduction; //對(duì)同一warp內(nèi)線程求得的中間結(jié)果進(jìn)行并行歸約求和

If (lane==0)y[row]+=vals[threadIdx.x]; //每個(gè)warp內(nèi)的第一個(gè)線程保存結(jié)果

end.

3.4 向量?jī)?nèi)積并行化

為了提高內(nèi)積計(jì)算的并行度, 本文使用一種CPU+GPU的混合歸約方法, 核心思想是歸約操作分兩步完成： 第一步將參與運(yùn)算的向量劃分成多個(gè)部分, 由每個(gè)block計(jì)算一對(duì)小向量的內(nèi)積, block中的每個(gè)線程負(fù)責(zé)對(duì)應(yīng)標(biāo)量的乘積, 計(jì)算結(jié)果存儲(chǔ)在共享內(nèi)存中, 然后將結(jié)果相加得到小向量的內(nèi)積； 第二步將第一步每個(gè)block求得的小向量?jī)?nèi)積傳送到CPU, 由CPU累加后得到最終結(jié)果.

由于block個(gè)數(shù)不能大于512, 為了充分利用GPU上的計(jì)算資源, 本文動(dòng)態(tài)分配每個(gè)block內(nèi)的線程數(shù)(即每個(gè)小向量的長(zhǎng)度). 小向量長(zhǎng)度m最小為128, 當(dāng)n>128×512時(shí),m=256； 當(dāng)n>256×512時(shí),m=512； 當(dāng)n>512×512時(shí),m=512×((n-1)/(512×512)+1). 這樣不同規(guī)模的問題都能得到最佳warp塊裝載量.

Block大小為128的向量?jī)?nèi)積Kernel.

thread_id=blockDim.x×blockIdx.x+threadIdx.x;； //每個(gè)線程的編號(hào)

tid=threadIdx.x//每個(gè)線程在block中的編號(hào)

_shared_res[128]; //定義共享內(nèi)存數(shù)組

res[tid]=v[thread_id]×w[thread_id]; //中間結(jié)果存入數(shù)組res中；

for (0≤tid<128)

Reduction； //對(duì)同一block內(nèi)線程求得的中間結(jié)果進(jìn)行并行歸約求和

if(tid<1)y[blockIdx.x]=res[0]+res[1]; //每個(gè)block內(nèi)的第一個(gè)線程保存結(jié)果

end.

3.5 向量加減、 向量乘標(biāo)量和預(yù)處理方程組求解

向量加減和向量乘標(biāo)量計(jì)算較簡(jiǎn)單, 為了減少CPU與GPU間的通信開銷, 向量加減和向量乘標(biāo)量計(jì)算也在GPU上完成. 本文借鑒CUDA SDK中的方法實(shí)現(xiàn). 對(duì)于求解預(yù)處理方程組時(shí)遇到的兩次三角方程組求解由CPU完成.

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

4.1 實(shí)驗(yàn)平臺(tái)

本文實(shí)驗(yàn)在PC機(jī)上完成, CPU為Intel Core Quad CPU Q8300@ 2.50 GHz, 內(nèi)存2 GRAM, 顯卡芯片型號(hào)為Geforce GTX460, 操作系統(tǒng)為32位的CentOS 5.5, 使用Intel Fortran編譯器和Nvdia的nvcc編譯器對(duì)混合代碼進(jìn)行編譯, CUDA工具包的版本是3.0.

4.2 實(shí)驗(yàn)設(shè)計(jì)與結(jié)果分析

由實(shí)際應(yīng)用產(chǎn)生的稀疏矩陣規(guī)模分別為： 150×150, 1 068×1 068, 2 583×2 583, 6 453×6 453.

圖5 GPU并行程序各部分執(zhí)行時(shí)間百分?jǐn)?shù)Fig.5 Time proportion of program execution on GPU

實(shí)驗(yàn)1為了驗(yàn)證CPU到GPU的數(shù)據(jù)傳輸對(duì)并行算法執(zhí)行效率的影響, 本文針對(duì)不同規(guī)模的數(shù)據(jù)集進(jìn)行了測(cè)試. 實(shí)驗(yàn)結(jié)果如圖5所示.

由圖5可見, 當(dāng)稀疏矩陣規(guī)模為150×150時(shí), 從CPU到GPU的數(shù)據(jù)傳輸占據(jù)了程序執(zhí)行的大部分時(shí)間, 是影響程序執(zhí)行速度的主要瓶頸.程序中的數(shù)據(jù)傳輸主要發(fā)生在初始化階段, 即將以CSR格式存儲(chǔ)的系數(shù)矩陣從內(nèi)存一次性拷貝到GPU顯存的過程.

實(shí)驗(yàn)2為了驗(yàn)證本文提出的算法對(duì)SPMV和向量?jī)?nèi)積算法的加速效果, 使用CPU(Q8300)和GPU(GTX460)各執(zhí)行100次操作, 進(jìn)行性能對(duì)比, 結(jié)果如圖6和圖7所示.

圖6 矩陣向量乘執(zhí)行100次耗時(shí)Fig.6 Time spent on 100 times SPMV

圖7 向量?jī)?nèi)積操作執(zhí)行100次耗時(shí)Fig.7 Time spent on 100 times inner product

由圖6和圖7可見, 當(dāng)矩陣規(guī)模只有150×150時(shí), 稀疏矩陣向量乘計(jì)算在CPU和GPU上的執(zhí)行時(shí)間基本持平, 而GPU上的向量?jī)?nèi)積計(jì)算比CPU串行執(zhí)行更耗時(shí). 這是因?yàn)楫?dāng)矩陣規(guī)模太小時(shí), GPU上同時(shí)運(yùn)行的線程數(shù)有限, 導(dǎo)致流處理器大量空置, 相當(dāng)于串行執(zhí)行, 且還要考慮其中的同步開銷. 而隨著矩陣規(guī)模的增大, GPU上計(jì)算資源占有率提高, 優(yōu)勢(shì)逐漸顯現(xiàn). 當(dāng)矩陣規(guī)模達(dá)到6 453×6 453時(shí), GPU上執(zhí)行單次SPMV計(jì)算相對(duì)CPU的加速比達(dá)到3.4倍.

實(shí)驗(yàn)3為了驗(yàn)證本文工作對(duì)整個(gè)地下多相流數(shù)值模擬的加速效果, 本文將并行算法整合進(jìn)TOUGHREACT, 對(duì)整體執(zhí)行時(shí)間進(jìn)行測(cè)試, 測(cè)試結(jié)果列于表2.

表2 CPU/GPU數(shù)值模擬執(zhí)行時(shí)間對(duì)比

由表2可見, 隨著矩陣規(guī)模的增大, 求解所需的迭代次數(shù)越來越多. 在雙精度條件下, CPU程序和GPU程序都能得到較好的收斂效果(矩陣規(guī)模為1 068×1 068時(shí)因無法得到收斂解所以未列出). 隨著迭代次數(shù)的增加, GPU程序進(jìn)行數(shù)值模擬節(jié)省的絕對(duì)執(zhí)行時(shí)間越多, GPU并行程序的執(zhí)行效果越好.

綜上所述, 本文利用CUDA編程模型特殊的通用計(jì)算模式, 對(duì)TOUGHREACT軟件中的大規(guī)模稀疏線性方程組求解器進(jìn)行了并行化改進(jìn), 實(shí)現(xiàn)了基于GPU求解線性方程組的預(yù)處理共軛梯度法. 隨著方程組的規(guī)模及稀疏矩陣維度的增多, 基于GPU并行算法的單次執(zhí)行時(shí)間比CPU串行算法有1.7～3.4倍的加速比.

[1] XU Tian-fu, Sonnenthal E, Spycher N, et al. TOURGHREACT: A Simulation Program for Non-isothermal Multiphase Reactive Geochemical Transport in Variably Saturated Geologic Media [J]. Computers & Geosciences, 2006, 32(2): 145-165.

[2] WANG Ming-liang, Klie H, Parashar M, et al. Solving Sparse Linear Systems on NVIDIA Tesla GPUs [C]//ICCS’09 Proceedings of the 9th International Conference on Computational Science. Berlin: Springer-Verlag, 2009: 864-873.

[3] ZHANG Ke-ni, WU Yu-shu, Bodvarsson G S. Parallel Computing Simulation of Fluid Flow in the Unsaturated Zone of Yucca Mountain, Nevada [J]. Journal of Contaminant Hydrology, 2003(62/63): 381-399.

[4] Huang J, Christ J A, Goltz M N. An Assembly Model for Simulation of Large-Scale Ground Water Flow and Transport [J]. Ground Water, 2008, 46(6): 882-892.

[5] DONG Yan-hui, LI Guo-min. A Parallel PCG Solver for MODFLOW [J]. Ground Water, 2009, 47(6): 845-850.

[6] JI Xiao-hui, CHENG Tang-pei, WANG Qun. CUDA-Based Solver for Large-Scale Groundwater Flow Simulation [J]. Engineering with Computers, 2012, 28(1): 13-19.

[7] Narasimhan T N, Witherspoon P A. An Integrated Finite Difference Method for Analyzing Fluid Flow in Porous Media [J]. Water Resources Research, 1976, 12(1): 57-64.

[8] Pruess K. TOUGH2: A General-Purpose Numerical Simulator for Multiphase Fluid and Heart Flow [R]. Berkeley: Lawrence Berkely Laboratory, 1991.

[9] ZHANG Ke-ni, WU Yu-shu, Ding C, et al. Parallel Computing Techniques for Large-Scale Reservoir Simulation of Multi-component and Multiphase Fluid Flow [C]//Proceedings of the 2001 SPE Reservoir Simulation Symposium. Houston, Texas: Society of Petroleum Engineers, 2001.

[10] Saad Y. Iterative Methods for Sparse Linear Systems [M]. Philadelphia: SIAM, 2003.

[11] Bell N, Garland M. Efficient Sparse Matrix Vector Multiplication on CUDA [EB]. 2008. http://mgarland.org.fields.papers.nvr-2008-004.pdf.