耦合溶質運移和化學反應模擬的并行優(yōu)化

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

(1. 吉林大學 計算機科學與技術學院,長春 130012;2. 吉林大學 地下能源及廢物處置研究所,長春 130021)

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

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

1 溶質運移化學耦合過程模型

1.1 基本過程

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

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

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

1.2 數學模型

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

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

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

2 耦合模擬過程性能分析

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

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

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

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

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

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

3 化學反應計算部分并行優(yōu)化

化學反應計算的復雜度直接與地質模型中的離散網格數相關,因此,需要對依次計算的順序進行拆分.

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

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

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

……

END DO

3.2 進程同步

化學反應的非線性方程組是對單一網格建立的,所以對于每個獨立進程,在化學反應部分不需要交換數據,可并行執(zhí)行. 但在溶質運移結束后、 化學反應計算前,要有一個同步通訊過程,原因如下：

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

3.3 并行化算法的實現

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

對耦合過程的并行化處理

溶質運移/化學反應CYCLE

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

Call Partition_job //劃分處理

Call Get_local_index //下標數組賦值

溶質運移……

地球化學反應模擬

對每個離散網格進行化學方程求解

Call Pack_array

Call Comm_intern //MPI同步操作

Call Unpack_array

DO (1,nprocs)

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

Call Assign //對網格施加參數信息

//濃度初始值處理

Call Newton_raphson //Newton法求解

…… //反應對水流狀態(tài)反饋處理

Call Update //對當前節(jié)點狀態(tài)進行更新

Call Conver //收斂檢驗處理

END DO

END DO

數據文件寫入

溶質運移/化學反應過程結束

4 實驗結果

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

4.1 算法執(zhí)行效率

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

圖4 50個網格的測試結果Fig.4 Results of input sample into 50 grids

圖5 861個網格的測試結果Fig.5 Results of input sample into 861 grids

4.2 算法的加速比

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

(2)

可知,上述3個實驗結果的加速比如圖7所示. 由圖7可見,化學反應的計算時間隨著分布在簡單機群上進程數的增加而縮短,隨著計算進程的增加,加速比的值逐漸上升,平均獲得了2～3倍的加速.

圖6 2 151個網格的測試結果Fig.6 Results of input sample into 2 151 grids

圖7 對不同測試用例的加速比Fig.7 Speed up on input samples

3個不同測試實例運行時間加速比較接近,在并行度相同的情況下,隨著網格數目的增加,加速比的值有所下降,這是因為化學狀態(tài)數組的維度直接與模擬網格的規(guī)模相關,拆分數組的聚合操作時間會相應地隨網格數目增加. 進而,對MPI操作緩沖區(qū)的數據量也有影響,導致MPI操作的開銷增加. 式(2)表明,并行開銷的增加會降低加速比. 在本文測試實例中,隨著問題規(guī)模的增大,并行開銷的增長速度略大于計算時間的增長,導致加速比曲線的斜率減小.

綜上可見,場地級耦合溶質運移與化學反應的應用需求,會極大增加數值模擬的計算時間. 需要在時間效率上對耦合數值模擬過程做出優(yōu)化. 本文通過分析耦合過程的數學模型及實際的執(zhí)行時間比例,針對占用較多計算時間的化學反應模擬部分設計了并行優(yōu)化方法. 并對不同規(guī)模的輸入模型在測試集群上進行了實驗. 實驗結果表明,該方法對耦合計算時間平均具有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.