魏晓辉,李维山,李洪亮,朱 彤,许天福

(1. 吉林大学 计算机科学与技术学院,长春 130012;2. 吉林大学 地下能源及废物处置研究所,长春 130021)

TOUGHREACT[1]作为在渗透性裂隙介质内非等温多相流的化学反应模拟程序,是在多相水流和热流模拟器TOUGH2[2-8]上加入化学反应及溶质运移的模拟过程而开发的,可用于模拟一维、 二维和三维孔隙或裂隙物理或化学环境复杂耦合过程的模拟. 根据模拟计算机平台的计算性能,可对不同种类、 不同数目的化学物质在液态、 气态及固态中反应过程进行数值模拟. TOUGHREACT在二氧化碳地址封存、 核废料处置、 地热能源开发及环境修复等领域应用广泛. 随着如场地级二氧化碳地质封存等应用[9]需求的不断增多,地质模型的复杂度和尺度也明显加大. 目前,串行版TOUGHREACT在处理这些新需求时显现了计算能力上的局限,需要提高计算效率[10-20].

本文通过分析耦合地球化学反应过程的特点及程序执行时间上的热点,针对分析结果设计并实现了并行求解算法框架. 通过实际地质模拟模型对程序改进后的执行效果进行了实验. 实验结果表明,该方法对耦合模拟计算具有良好的加速效果,在测试机群上可以有2～3倍的加速比.

1 溶质运移化学耦合过程模型

1.1 基本过程

TOUGHREACT在离散方法上采用与TOUGH2相同的积分有限差分法(integral finite difference method,IFDM),对流体和热量运移过程的模拟与TOUGH2一致,并加入水溶液和气态的溶质运移及化学反应过程模拟.

图1 TOUGHREACT基本过程模拟Fig.1 Main procedure of simulation via TOUGHREACT

TOUGHREACT采用不同过程间顺次数值计算处理物理化学耦合的数值模拟过程：在完成对水流方程的求解后,水流速度与相饱和度的数值被用于溶质运移的模拟过程. 对溶质运移的数值过程求解则按各组分顺序进行. 在耦合过程中,从求解溶质运移方程得到的结果浓度值被代入化学反应模型中. 该混合动力学平衡的化学反应系统方程是在网格块顺序下,依次用Newton-Raphson方法求解. 该化学和溶质运移的过程被迭代进行,直到耦合系统收敛. 如图1所示,虚框内为溶质运移和化学反应式耦合模拟部分.

1.2 数学模型

TOUGHREACT的耦合运移化学反应数学模型[21]以守恒方程形式给出,对大多数化学组分,溶质运移在液态中进行. 在顺序迭代法(SIA)中,质量运移方程组和化学反应方程组被视为两个相对独立的子系统,在求解过程中以一种顺序迭代的次序分别求解. 当在液态状态下的化学反应被认定为局部平衡状态后,质量运移方程即可用总体溶解组分浓度的形式表示. 通过对各种质量积分项的求和,可得在液态中多组分的化学溶质运移方程组:

对气态组分浓度,方程的基本形式一样,区别仅在于浓度表达式的不同,所以求解方法与液态运移方程组相同.

混合动力学平衡化学系统的方程组建是基于基本组分质量守恒和化学平衡. 在数值解法上,与溶质运移的求解方法相同,对每个网格上建立的非线性方程也采用Newton-Raphson方法求解. 通过软件中参数的预设值,判断该迭代过程是否收敛. 化学反应计算部分的过程对每个网格相对独立,计算时不需交换信息,所以,本文在该计算部分进行并行优化.

2 耦合模拟过程性能分析

图2 计算过程的循环Fig.2 Cycle of simulation procedure

在软件设计上TOUGHREACT基本沿用了TOUGH2的架构,时间步的控制是最大的循环体,循环体内部是主要的模拟过程. 多相水流模拟和耦合溶质运移化学反应是相对独立的过程,可视为依次执行,每个部分都有独立的收敛性判断及状态变量更新. TOUGHREACT在耦合计算部分的迭代流程如图2所示.

从规模上看,溶质运移和化学反应两个方程系统的复杂度基本相同,但形式不同. 溶质运移一次建立针对所有网格的方程;化学反应是每次针对单一网格组建方程,但有网格规模大小的循环复杂度. 为了解程序在运行时的热点,下面用不同规模模拟问题作为输入实例,对程序执行情况进行分析.

图3 各部分执行时间比例Fig.3 Executing time proportion

一般各计算部分的时间比例随着输入数据的不同而变化,例如对地质模型的分层和网格属性数据值的变化,会影响到具体执行时间. 本文预先测试了各计算部分所占用的时间,分别用不同规模(861,2 151,7 550)的离散地质模型作为输入,对原有串行程序的执行时间做出分析,如图3所示. 运行的结果数据显示,多相水流和化学反应的时间比例较大,占执行时间的主要部分,而溶质运移的时间比例较小. 溶质运移与化学反应所占时间的比值,随着模拟问题规模的不同有所变化,但占用时间单位的量级比约为1∶10.

在耦合化学反应过程中,化学反应模拟计算是计算任务最密集的部分. 对化学反应过程的并行化,在执行并行度小于10的情况下,可有效缩短总体耦合过程的执行时间. 因此本文对耦合溶质运移和化学反应的并行化优化工作,主要针对化学反应的并行化.

3 化学反应计算部分并行优化

化学反应计算的复杂度直接与地质模型中的离散网格数相关,因此,需要对依次计算的顺序进行拆分.

3.1 循环分割计算

对于化学计算的循环部分,假设需要模拟的离散网格数目为NNEL,并行进程数为nprocs,平均每个进程负责计算的局部网格数目为NNEL/nprocs,进程负载最多网格数目为NNEL/nprocs+mod(NNEL,nprocs). 每个进程在执行代码上复用对单个网格单元的化学方程系数组建,Newton-Raphson迭代求解,化学组分状态变量更新等过程代码. 在子任务的循环下标处理方面,可用局部下标数组实现,如果进程的标识变量为myid,则代码可按如下方式实现：

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

……

END DO

3.2 进程同步

化学反应的非线性方程组是对单一网格建立的,所以对于每个独立进程,在化学反应部分不需要交换数据,可并行执行. 但在溶质运移结束后、 化学反应计算前,要有一个同步通讯过程,原因如下：

1) 对化学部分计算的拆分,不能保证划分的任务绝对均衡,且受各网格方程子系统Newton-Raphson迭代收敛速度的影响,各子进程的计算完成时间也会不同;

3.3 并行化算法的实现

为了能在分布式内存集群上进行并行化计算,本文采用MPI库实现不同计算节点间的通讯. 为满足分析的要求,在耦合过程的并行化实现了如下过程：通讯环境初始化的Comm_init子例程(封装了MPI环境初始化操作);对网格计算任务进行划分处理的Partition_job子例程;对数组的打包处理Pack_array和解包Unpack_array子例程;负责同步通讯的Comm_intern子例程(封装了MPI的通讯操作)及对化学循环的局部下标Local_lower和上标Local_upper处理的Get_local_index等子例程. 并行化后的程序描述如下:

对耦合过程的并行化处理

溶质运移/化学反应CYCLE

Call Comm_init //MPI环境初始化

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法求解

…… //反应对水流状态反馈处理

Call Update //对当前节点状态进行更新

Call Conver //收敛检验处理

END DO

END DO

数据文件写入

溶质运移/化学反应过程结束

4 实验结果

实验测试平台由4台PC服务器组成(CPU Intel Core2 Duo E7400 2.80 GHz,1 Gb内存)；服务器之间采用百兆以太网连接. 测试环境: Linux操作系统(Cent OS 5.5);MPI并行运行环境库(MPICH2 1.3.1);Fortran编译器(Intel Fortran Compiler 11.1 Linux).

4.1 算法执行效率

为验证本文提出的并行优化方法能有效加速耦合过程计算,使用3个带耦合化学反应过程地质模型实例. 在输入文件中MESH文件定义的网格规模分别为50,861,2 151,测试结果分别如图4～图6所示. 由图4～图6可见,对化学计算部分的并行化处理可有效加速耦合化学计算过程. 溶质运移的计算部分在程序改进后,对每个进程使用相同的代码和数据,所以该部分的计算时间较稳定.

图4 50个网格的测试结果Fig.4 Results of input sample into 50 grids

图5 861个网格的测试结果Fig.5 Results of input sample into 861 grids

4.2 算法的加速比

若串行执行时间为Ts,并行执行时间为Tp,并行执行进程数为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个不同测试实例运行时间加速比较接近,在并行度相同的情况下,随着网格数目的增加,加速比的值有所下降,这是因为化学状态数组的维度直接与模拟网格的规模相关,拆分数组的聚合操作时间会相应地随网格数目增加. 进而,对MPI操作缓冲区的数据量也有影响,导致MPI操作的开销增加. 式(2)表明,并行开销的增加会降低加速比. 在本文测试实例中,随着问题规模的增大,并行开销的增长速度略大于计算时间的增长,导致加速比曲线的斜率减小.

综上可见,场地级耦合溶质运移与化学反应的应用需求,会极大增加数值模拟的计算时间. 需要在时间效率上对耦合数值模拟过程做出优化. 本文通过分析耦合过程的数学模型及实际的执行时间比例,针对占用较多计算时间的化学反应模拟部分设计了并行优化方法. 并对不同规模的输入模型在测试集群上进行了实验. 实验结果表明,该方法对耦合计算时间平均具有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.