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

(1. 吉林大學(xué) 計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院,長(zhǎng)春 130012;2. 吉林大學(xué) 地下能源及廢物處置研究所,長(zhǎng)春 130021)

TOUGHREACT[1]作為在滲透性裂隙介質(zhì)內(nèi)非等溫多相流的化學(xué)反應(yīng)模擬程序,是在多相水流和熱流模擬器TOUGH2[2-8]上加入化學(xué)反應(yīng)及溶質(zhì)運(yùn)移的模擬過程而開發(fā)的,可用于模擬一維、 二維和三維孔隙或裂隙物理或化學(xué)環(huán)境復(fù)雜耦合過程的模擬. 根據(jù)模擬計(jì)算機(jī)平臺(tái)的計(jì)算性能,可對(duì)不同種類、 不同數(shù)目的化學(xué)物質(zhì)在液態(tài)、 氣態(tài)及固態(tài)中反應(yīng)過程進(jìn)行數(shù)值模擬. TOUGHREACT在二氧化碳地址封存、 核廢料處置、 地?zé)崮茉撮_發(fā)及環(huán)境修復(fù)等領(lǐng)域應(yīng)用廣泛. 隨著如場(chǎng)地級(jí)二氧化碳地質(zhì)封存等應(yīng)用[9]需求的不斷增多,地質(zhì)模型的復(fù)雜度和尺度也明顯加大. 目前,串行版TOUGHREACT在處理這些新需求時(shí)顯現(xiàn)了計(jì)算能力上的局限,需要提高計(jì)算效率[10-20].

本文通過分析耦合地球化學(xué)反應(yīng)過程的特點(diǎn)及程序執(zhí)行時(shí)間上的熱點(diǎn),針對(duì)分析結(jié)果設(shè)計(jì)并實(shí)現(xiàn)了并行求解算法框架. 通過實(shí)際地質(zhì)模擬模型對(duì)程序改進(jìn)后的執(zhí)行效果進(jìn)行了實(shí)驗(yàn). 實(shí)驗(yàn)結(jié)果表明,該方法對(duì)耦合模擬計(jì)算具有良好的加速效果,在測(cè)試機(jī)群上可以有2～3倍的加速比.

1 溶質(zhì)運(yùn)移化學(xué)耦合過程模型

1.1 基本過程

TOUGHREACT在離散方法上采用與TOUGH2相同的積分有限差分法(integral finite difference method,IFDM),對(duì)流體和熱量運(yùn)移過程的模擬與TOUGH2一致,并加入水溶液和氣態(tài)的溶質(zhì)運(yùn)移及化學(xué)反應(yīng)過程模擬.

圖1 TOUGHREACT基本過程模擬Fig.1 Main procedure of simulation via TOUGHREACT

TOUGHREACT采用不同過程間順次數(shù)值計(jì)算處理物理化學(xué)耦合的數(shù)值模擬過程：在完成對(duì)水流方程的求解后,水流速度與相飽和度的數(shù)值被用于溶質(zhì)運(yùn)移的模擬過程. 對(duì)溶質(zhì)運(yùn)移的數(shù)值過程求解則按各組分順序進(jìn)行. 在耦合過程中,從求解溶質(zhì)運(yùn)移方程得到的結(jié)果濃度值被代入化學(xué)反應(yīng)模型中. 該混合動(dòng)力學(xué)平衡的化學(xué)反應(yīng)系統(tǒng)方程是在網(wǎng)格塊順序下,依次用Newton-Raphson方法求解. 該化學(xué)和溶質(zhì)運(yùn)移的過程被迭代進(jìn)行,直到耦合系統(tǒng)收斂. 如圖1所示,虛框內(nèi)為溶質(zhì)運(yùn)移和化學(xué)反應(yīng)式耦合模擬部分.

1.2 數(shù)學(xué)模型

TOUGHREACT的耦合運(yùn)移化學(xué)反應(yīng)數(shù)學(xué)模型[21]以守恒方程形式給出,對(duì)大多數(shù)化學(xué)組分,溶質(zhì)運(yùn)移在液態(tài)中進(jìn)行. 在順序迭代法(SIA)中,質(zhì)量運(yùn)移方程組和化學(xué)反應(yīng)方程組被視為兩個(gè)相對(duì)獨(dú)立的子系統(tǒng),在求解過程中以一種順序迭代的次序分別求解. 當(dāng)在液態(tài)狀態(tài)下的化學(xué)反應(yīng)被認(rèn)定為局部平衡狀態(tài)后,質(zhì)量運(yùn)移方程即可用總體溶解組分濃度的形式表示. 通過對(duì)各種質(zhì)量積分項(xiàng)的求和,可得在液態(tài)中多組分的化學(xué)溶質(zhì)運(yùn)移方程組:

對(duì)氣態(tài)組分濃度,方程的基本形式一樣,區(qū)別僅在于濃度表達(dá)式的不同,所以求解方法與液態(tài)運(yùn)移方程組相同.

混合動(dòng)力學(xué)平衡化學(xué)系統(tǒng)的方程組建是基于基本組分質(zhì)量守恒和化學(xué)平衡. 在數(shù)值解法上,與溶質(zhì)運(yùn)移的求解方法相同,對(duì)每個(gè)網(wǎng)格上建立的非線性方程也采用Newton-Raphson方法求解. 通過軟件中參數(shù)的預(yù)設(shè)值,判斷該迭代過程是否收斂. 化學(xué)反應(yīng)計(jì)算部分的過程對(duì)每個(gè)網(wǎng)格相對(duì)獨(dú)立,計(jì)算時(shí)不需交換信息,所以,本文在該計(jì)算部分進(jìn)行并行優(yōu)化.

2 耦合模擬過程性能分析

圖2 計(jì)算過程的循環(huán)Fig.2 Cycle of simulation procedure

在軟件設(shè)計(jì)上TOUGHREACT基本沿用了TOUGH2的架構(gòu),時(shí)間步的控制是最大的循環(huán)體,循環(huán)體內(nèi)部是主要的模擬過程. 多相水流模擬和耦合溶質(zhì)運(yùn)移化學(xué)反應(yīng)是相對(duì)獨(dú)立的過程,可視為依次執(zhí)行,每個(gè)部分都有獨(dú)立的收斂性判斷及狀態(tài)變量更新. TOUGHREACT在耦合計(jì)算部分的迭代流程如圖2所示.

從規(guī)模上看,溶質(zhì)運(yùn)移和化學(xué)反應(yīng)兩個(gè)方程系統(tǒng)的復(fù)雜度基本相同,但形式不同. 溶質(zhì)運(yùn)移一次建立針對(duì)所有網(wǎng)格的方程;化學(xué)反應(yīng)是每次針對(duì)單一網(wǎng)格組建方程,但有網(wǎng)格規(guī)模大小的循環(huán)復(fù)雜度. 為了解程序在運(yùn)行時(shí)的熱點(diǎn),下面用不同規(guī)模模擬問題作為輸入實(shí)例,對(duì)程序執(zhí)行情況進(jìn)行分析.

圖3 各部分執(zhí)行時(shí)間比例Fig.3 Executing time proportion

一般各計(jì)算部分的時(shí)間比例隨著輸入數(shù)據(jù)的不同而變化,例如對(duì)地質(zhì)模型的分層和網(wǎng)格屬性數(shù)據(jù)值的變化,會(huì)影響到具體執(zhí)行時(shí)間. 本文預(yù)先測(cè)試了各計(jì)算部分所占用的時(shí)間,分別用不同規(guī)模(861,2 151,7 550)的離散地質(zhì)模型作為輸入,對(duì)原有串行程序的執(zhí)行時(shí)間做出分析,如圖3所示. 運(yùn)行的結(jié)果數(shù)據(jù)顯示,多相水流和化學(xué)反應(yīng)的時(shí)間比例較大,占執(zhí)行時(shí)間的主要部分,而溶質(zhì)運(yùn)移的時(shí)間比例較小. 溶質(zhì)運(yùn)移與化學(xué)反應(yīng)所占時(shí)間的比值,隨著模擬問題規(guī)模的不同有所變化,但占用時(shí)間單位的量級(jí)比約為1∶10.

在耦合化學(xué)反應(yīng)過程中,化學(xué)反應(yīng)模擬計(jì)算是計(jì)算任務(wù)最密集的部分. 對(duì)化學(xué)反應(yīng)過程的并行化,在執(zhí)行并行度小于10的情況下,可有效縮短總體耦合過程的執(zhí)行時(shí)間. 因此本文對(duì)耦合溶質(zhì)運(yùn)移和化學(xué)反應(yīng)的并行化優(yōu)化工作,主要針對(duì)化學(xué)反應(yīng)的并行化.

3 化學(xué)反應(yīng)計(jì)算部分并行優(yōu)化

化學(xué)反應(yīng)計(jì)算的復(fù)雜度直接與地質(zhì)模型中的離散網(wǎng)格數(shù)相關(guān),因此,需要對(duì)依次計(jì)算的順序進(jìn)行拆分.

3.1 循環(huán)分割計(jì)算

對(duì)于化學(xué)計(jì)算的循環(huán)部分,假設(shè)需要模擬的離散網(wǎng)格數(shù)目為NNEL,并行進(jìn)程數(shù)為nprocs,平均每個(gè)進(jìn)程負(fù)責(zé)計(jì)算的局部網(wǎng)格數(shù)目為NNEL/nprocs,進(jìn)程負(fù)載最多網(wǎng)格數(shù)目為NNEL/nprocs+mod(NNEL,nprocs). 每個(gè)進(jìn)程在執(zhí)行代碼上復(fù)用對(duì)單個(gè)網(wǎng)格單元的化學(xué)方程系數(shù)組建,Newton-Raphson迭代求解,化學(xué)組分狀態(tài)變量更新等過程代碼. 在子任務(wù)的循環(huán)下標(biāo)處理方面,可用局部下標(biāo)數(shù)組實(shí)現(xiàn),如果進(jìn)程的標(biāo)識(shí)變量為myid,則代碼可按如下方式實(shí)現(xiàn)：

Do (Local_lower[myid],Local_upper[myid])

……

END DO

3.2 進(jìn)程同步

化學(xué)反應(yīng)的非線性方程組是對(duì)單一網(wǎng)格建立的,所以對(duì)于每個(gè)獨(dú)立進(jìn)程,在化學(xué)反應(yīng)部分不需要交換數(shù)據(jù),可并行執(zhí)行. 但在溶質(zhì)運(yùn)移結(jié)束后、 化學(xué)反應(yīng)計(jì)算前,要有一個(gè)同步通訊過程,原因如下：

1) 對(duì)化學(xué)部分計(jì)算的拆分,不能保證劃分的任務(wù)絕對(duì)均衡,且受各網(wǎng)格方程子系統(tǒng)Newton-Raphson迭代收斂速度的影響,各子進(jìn)程的計(jì)算完成時(shí)間也會(huì)不同;

3.3 并行化算法的實(shí)現(xiàn)

為了能在分布式內(nèi)存集群上進(jìn)行并行化計(jì)算,本文采用MPI庫(kù)實(shí)現(xiàn)不同計(jì)算節(jié)點(diǎn)間的通訊. 為滿足分析的要求,在耦合過程的并行化實(shí)現(xiàn)了如下過程：通訊環(huán)境初始化的Comm_init子例程(封裝了MPI環(huán)境初始化操作);對(duì)網(wǎng)格計(jì)算任務(wù)進(jìn)行劃分處理的Partition_job子例程;對(duì)數(shù)組的打包處理Pack_array和解包Unpack_array子例程;負(fù)責(zé)同步通訊的Comm_intern子例程(封裝了MPI的通訊操作)及對(duì)化學(xué)循環(huán)的局部下標(biāo)Local_lower和上標(biāo)Local_upper處理的Get_local_index等子例程. 并行化后的程序描述如下:

對(duì)耦合過程的并行化處理

溶質(zhì)運(yùn)移/化學(xué)反應(yīng)CYCLE

Call Comm_init //MPI環(huán)境初始化

Call Partition_job //劃分處理

Call Get_local_index //下標(biāo)數(shù)組賦值

溶質(zhì)運(yùn)移……

地球化學(xué)反應(yīng)模擬

對(duì)每個(gè)離散網(wǎng)格進(jìn)行化學(xué)方程求解

Call Pack_array

Call Comm_intern //MPI同步操作

Call Unpack_array

DO (1,nprocs)

DO (Local_lower[myid],Local_upper[myid])

Call Assign //對(duì)網(wǎng)格施加參數(shù)信息

//濃度初始值處理

Call Newton_raphson //Newton法求解

…… //反應(yīng)對(duì)水流狀態(tài)反饋處理

Call Update //對(duì)當(dāng)前節(jié)點(diǎn)狀態(tài)進(jìn)行更新

Call Conver //收斂檢驗(yàn)處理

END DO

END DO

數(shù)據(jù)文件寫入

溶質(zhì)運(yùn)移/化學(xué)反應(yīng)過程結(jié)束

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

實(shí)驗(yàn)測(cè)試平臺(tái)由4臺(tái)PC服務(wù)器組成(CPU Intel Core2 Duo E7400 2.80 GHz,1 Gb內(nèi)存)；服務(wù)器之間采用百兆以太網(wǎng)連接. 測(cè)試環(huán)境: Linux操作系統(tǒng)(Cent OS 5.5);MPI并行運(yùn)行環(huán)境庫(kù)(MPICH2 1.3.1);Fortran編譯器(Intel Fortran Compiler 11.1 Linux).

4.1 算法執(zhí)行效率

為驗(yàn)證本文提出的并行優(yōu)化方法能有效加速耦合過程計(jì)算,使用3個(gè)帶耦合化學(xué)反應(yīng)過程地質(zhì)模型實(shí)例. 在輸入文件中MESH文件定義的網(wǎng)格規(guī)模分別為50,861,2 151,測(cè)試結(jié)果分別如圖4～圖6所示. 由圖4～圖6可見,對(duì)化學(xué)計(jì)算部分的并行化處理可有效加速耦合化學(xué)計(jì)算過程. 溶質(zhì)運(yùn)移的計(jì)算部分在程序改進(jìn)后,對(duì)每個(gè)進(jìn)程使用相同的代碼和數(shù)據(jù),所以該部分的計(jì)算時(shí)間較穩(wěn)定.

圖4 50個(gè)網(wǎng)格的測(cè)試結(jié)果Fig.4 Results of input sample into 50 grids

圖5 861個(gè)網(wǎng)格的測(cè)試結(jié)果Fig.5 Results of input sample into 861 grids

4.2 算法的加速比

若串行執(zhí)行時(shí)間為Ts,并行執(zhí)行時(shí)間為Tp,并行執(zhí)行進(jìn)程數(shù)為n,并行開銷為To,根據(jù)加速比的定義

(2)

可知,上述3個(gè)實(shí)驗(yàn)結(jié)果的加速比如圖7所示. 由圖7可見,化學(xué)反應(yīng)的計(jì)算時(shí)間隨著分布在簡(jiǎn)單機(jī)群上進(jìn)程數(shù)的增加而縮短,隨著計(jì)算進(jìn)程的增加,加速比的值逐漸上升,平均獲得了2～3倍的加速.

圖6 2 151個(gè)網(wǎng)格的測(cè)試結(jié)果Fig.6 Results of input sample into 2 151 grids

圖7 對(duì)不同測(cè)試用例的加速比Fig.7 Speed up on input samples

3個(gè)不同測(cè)試實(shí)例運(yùn)行時(shí)間加速比較接近,在并行度相同的情況下,隨著網(wǎng)格數(shù)目的增加,加速比的值有所下降,這是因?yàn)榛瘜W(xué)狀態(tài)數(shù)組的維度直接與模擬網(wǎng)格的規(guī)模相關(guān),拆分?jǐn)?shù)組的聚合操作時(shí)間會(huì)相應(yīng)地隨網(wǎng)格數(shù)目增加. 進(jìn)而,對(duì)MPI操作緩沖區(qū)的數(shù)據(jù)量也有影響,導(dǎo)致MPI操作的開銷增加. 式(2)表明,并行開銷的增加會(huì)降低加速比. 在本文測(cè)試實(shí)例中,隨著問題規(guī)模的增大,并行開銷的增長(zhǎng)速度略大于計(jì)算時(shí)間的增長(zhǎng),導(dǎo)致加速比曲線的斜率減小.

綜上可見,場(chǎng)地級(jí)耦合溶質(zhì)運(yùn)移與化學(xué)反應(yīng)的應(yīng)用需求,會(huì)極大增加數(shù)值模擬的計(jì)算時(shí)間. 需要在時(shí)間效率上對(duì)耦合數(shù)值模擬過程做出優(yōu)化. 本文通過分析耦合過程的數(shù)學(xué)模型及實(shí)際的執(zhí)行時(shí)間比例,針對(duì)占用較多計(jì)算時(shí)間的化學(xué)反應(yīng)模擬部分設(shè)計(jì)了并行優(yōu)化方法. 并對(duì)不同規(guī)模的輸入模型在測(cè)試集群上進(jìn)行了實(shí)驗(yàn). 實(shí)驗(yàn)結(jié)果表明,該方法對(duì)耦合計(jì)算時(shí)間平均具有2～3倍的加速效果.

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

[2] Pruess K. TOUGH 2: A General-Purpose Numerical Simulator for Multiphase Fluid and Heart Flow [R]. Berkeley: CA,1991.

[3] Elmroth E,Ding C,WU Yu-shu,et al. High Performance Computations for Large-Scale Simulations of Subsurface Multiphase Fluid and Heat Flow [J]. J Supercomputing,2001,18(3): 233-258.

[4] 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: SPE,2001: SPE66343.

[5] ZHANG Ke-ni,WU Yu-shu,Ding C,et al. TOUGH2_MP: A Parallel Version of TOUGH2 [C]//Proceedings of TOUGH Symposium 2003. Berkeley,California: Lawrence Berkeley National Laboratory,2003: 12-14.

[6] WU Yu-shu,ZHANG Ke-ni,Ding C,et al. An Effcient Parallel-Computing Method for Modeling Nonisothermal Multiphase Flow and Multicomponent Transport in Porous and Fractured Media [J]. Advances in Water Resources,2002,25(3): 243-261.

[7] ZHANG Ke-ni,Moridis G J,WU Yu-shu,et al. A Domain Decomposition Approach for Large Scale Simulations of Flow Processes in Hydrate Bearing Geologic Media [C/OL]. 2009-03-11. http://www.escholarshiporg/uc/item/4595r17h.

[8] Bhogeswara R,Killough J E. Parallel Linear Solvers for Reservoir Simulation: Generic Approach for Existing and Emerging Computer Architectures [J]. SPE Computer Applications,1994,6(1): 5-11.

[9] Audigane P,Gaus I,Czernichowski-Lauriol I,et al. Two-Dimensional Reactive Transport Modeling of CO2Injection in a Saline Aquifer at the Sleipner Site [J]. American Journal of Science,2007,307(7): 974-1008.

[10] Briens F J L,Wu C H,Gazdag J,et al. Compositional Reservoir Simulation in Parallel Supercomputing Environments [C]//Proceedings of 11th SPE Symposium on Reservoir Simulation. Anaheim,CA: SPE,1991: 125-133.

[11] Barua J,Horne R N. Improving the Performance of Parallel (and Series) Reservoir Simulation [C]//Proceedings of 10th SPE Symposium on Reservoir Simulation. Houston,TX: SPE,1989: 7-18.

[12] Meijerink J A,Daalen D T,Van,Hoogerbrugge P J,et al. Towards a More Efficient Parallel Reservoir Simulator [C]//Proceedings of 11th SPE Symposium on Reservoir Simulation. Anaheim,CA: SPE,1991: 107-116.

[13] Quandalle P,Moriano S. Vectorization and Parallel Processing of Models with Local Refinement [J]. SPE Advanced Technology Series,1993,1(2): 93-99.

[14] Wallis J R,Foster J A,Kendall R P. A New Parallel Iterative Linearsolution Method for Large-Scale Reservoir Simulation [C]//Proceedings of 11th SPE Symposium on Reservoir Simulation. Anaheim,CA: SPE,1991: 83-92.

[15] Chien M C H,Northrup E J. Vectorization and Parallel Processing of Local Grid Refinement and Adaptive Implicit Schemes in a General Purpose Reservoir Simulator [C]//Proceedings of 12th SPE Symposium on Reservoir Simulation. New Orleans,LA: SPE,1993: 279-290.

[16] Killough J E,Bhogeswara R. Simulation of Compositional Reservoir Phenomena on a Distributed Memory Parallel Computer [J]. Journal of Petroleum Technoligy,1991,43(11): 1368-1374.

[17] Wheeler J A,Smith R A. Reservoir Simulation on a Hypercube [J]. SPE Reservoir Eng,1990,5(4): 544-548.

[18] Kohar G,Killough J E. An Asynchronous Parallel Linear Equation Solution Technique [C]//Proceedings of 13th SPE Symposium on Reservoir Simulation. San Antonio,TX: SPE,1995: 507-520.

[19] Wang P,Balay S,Sepehrnoori K,et al. A Fully Implicit Parallel EOS Compositional Simulator for Large Scale Reservoir Simulation [C]//Proceedings of 15th SPE Symposium on Reservoir Simulation. Houston,TX: SPE,1999: 63-71.

[20] Vertiere S,Quettier L,Samier P,et al. Application of a Parallel Simulator to Industrial Test Cases [C]//Proceedings of 15th SPE Symposium on Reservoir Simulation. Houston,TX: SPE,1999: 93-105.

[21] XU Tian-fu,Sonnenthal E,Spycher N,et al. TOUGHREACT User’s Guide: A Simulation Program for Non-isothermal Multiphase Reactive Geochemical Transport in Variably Saturated Geologic Media [DB/OL]. 2004-05-24. http://escholarship.org/uc/item/8d43d056.