OpenMP在CO2地質(zhì)儲(chǔ)存數(shù)值模擬并行計(jì)算中的應(yīng)用

楊艷林， 許天福， 靖 晶

(1.武漢地質(zhì)調(diào)查中心，湖北 武漢 430205；2.吉林大學(xué)環(huán)境與資源學(xué)院，吉林 長(zhǎng)春 130021；3.中國(guó)地質(zhì)大學(xué)(武漢)環(huán)境學(xué)院, 湖北 武漢 430074)

CO2地質(zhì)封存作為CO2減排的重要方式，受到世界各國(guó)政府與大量學(xué)者的廣泛關(guān)注。數(shù)值模型在CO2地質(zhì)封存場(chǎng)地選擇、封存能力與安全性評(píng)估等方面發(fā)揮了重要作用[1～6]。TOUGH2-ECO2N程序是專門用于CO2地質(zhì)儲(chǔ)存的數(shù)值模擬器，模擬結(jié)果準(zhǔn)確可靠，得到了廣泛應(yīng)用[7～9]；但隨著研究的不斷深入，其對(duì)場(chǎng)地級(jí)精細(xì)模擬和長(zhǎng)時(shí)間跨度的高仿真數(shù)值模擬提出了迫切需求?；诖?，國(guó)內(nèi)外學(xué)者提出了OpenMP、PVM、MPI和GPU等并行計(jì)算技術(shù)[10]，同時(shí)開發(fā)出一些通用的并行求解器，如AZTEC、PETSc、Hypre[11]；以及區(qū)域分解工具，如Metis、ParMETIS。在CO2地質(zhì)封存領(lǐng)域也取得了一些成果，如國(guó)外PFLOTRAN、AM、HBGC123D以及PARSWMS等并行模擬器，成功應(yīng)用于CO2地質(zhì)封存；以及國(guó)內(nèi)張可霓等[12]基于MPI開發(fā)了TOUGH2-MP，可快速地完成CO2地質(zhì)封存數(shù)值模擬；朱彤[13]利用GPU，米東等[14]利用MPI加速TOUGHREACT，完成了注入CO2的水巖相互作用模擬。這些研究在加快計(jì)算速度和提高求解規(guī)模方面取得了一定的進(jìn)展，同時(shí)也為本文的研究提供了參考。

OpenMP是共享編程的典型編程模型，是計(jì)算機(jī)硬件、軟件廠商聯(lián)合發(fā)表的共享內(nèi)存編程應(yīng)用程序接口，具有編程簡(jiǎn)單、靈活、開發(fā)周期短、并行效率高、支持多種語言和多種操作系統(tǒng)等特點(diǎn)。本文在Windows PC機(jī)上，提出了利用OpenMP并行化TOUGH2-ECO2N模擬器的方法，大大增強(qiáng)了其計(jì)算規(guī)模及效率的能力。

1 TOUGH2-ECO2N數(shù)值模擬器簡(jiǎn)介

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

深部巖土體介質(zhì)中地下流體的控制方程是依據(jù)質(zhì)量和能量守恒準(zhǔn)則建立的，如下：

(1)

其中，Aκ可分為質(zhì)量項(xiàng)與能量項(xiàng)，質(zhì)量項(xiàng)主要包括水、鹽和CO2等，即式(2)；能量項(xiàng)是固相、液相與氣相中的能量和，即式(3)。同樣Fκ也分為質(zhì)量流量與能量流量，質(zhì)量流量主要采用達(dá)西定律來描述，即式(4)；能量流量采用熱傳導(dǎo)定律來刻畫，即式(5)。式中各參數(shù)意義見表1。

(2)

(3)

(4)

(5)

表1 模型中所涉及的符號(hào)及意義Table 1 Symbols and meaning in these models

1.2 模擬器程序結(jié)構(gòu)

TOUGH2是一個(gè)面向過程的程序，主體采用FORTRAN77語言編寫，包括一個(gè)主程序和若干個(gè)狀態(tài)方程模塊(圖1中的EOS，用于計(jì)算單元的次級(jí)主要變量)，不同狀態(tài)方程模塊，對(duì)應(yīng)不同功能的模擬器，如對(duì)于CO2地質(zhì)封存則采用其中的ECO2N模塊。通過對(duì)TOUGH2源程序分析，其程序流程見圖1。

圖1 TOUGH2模擬器程序簡(jiǎn)化流程圖Fig.1 Simplified flow chart of the TOUGH2 simulator

1.3 空間時(shí)間離散及矩陣存儲(chǔ)結(jié)構(gòu)

TOUGH2模擬器是采用積分有限差進(jìn)行空間離散、時(shí)間向前差分進(jìn)行時(shí)間離散，并通過程序中的“MULTI”來組建大規(guī)模矩陣方程(又稱Jacobian矩陣)。矩陣方程具有非對(duì)稱、強(qiáng)非線性和高稀疏的特點(diǎn)；最后采用三元數(shù)組的方式進(jìn)行存儲(chǔ)，生成后的矩陣形式見圖2。

圖2 Jacobian矩陣存儲(chǔ)(30個(gè)單元等溫模型)Fig.2 Jacobian matrix storage (Isothermal models of 30 cells)

1.4 模擬器耗時(shí)評(píng)估

為了了解程序運(yùn)行時(shí)的耗時(shí)部分，本文以中大規(guī)模的模型(2 500，10 000，40 000)為算例進(jìn)行分析。結(jié)合前面對(duì)TOUGH2源程序結(jié)構(gòu)的分析，可將模擬器分為四大過程：首先是對(duì)各種輸入文件進(jìn)行讀入；第二步是計(jì)算每個(gè)單元的狀態(tài)(EOS模塊)；第三步是組建矩陣方程(MULTI模塊)；最后是對(duì)矩陣方程進(jìn)行求解。第一部分耗時(shí)較小，第二、三和四部分是模擬器的主要耗時(shí)部分，其中方程求解所占的時(shí)間比例最大，達(dá)到70%以上，而且隨著模型規(guī)模的增大，這個(gè)比例還會(huì)增大(圖3)。

圖3 主要模塊執(zhí)行時(shí)間對(duì)比(20 a)Fig.3 Comparison of the module execution time (20 years)

單元的狀態(tài)更新與組建Jacobian矩陣所耗的時(shí)間相當(dāng)，因此在方程求解模塊進(jìn)行并行時(shí)需引起高度的重視。方程的求解方法，通常有直接法和迭代法，但由于直接法需要更多的計(jì)算機(jī)時(shí)和存儲(chǔ)空間，而迭代法彌補(bǔ)了其不足，尤其是處理大規(guī)模的矩陣時(shí)，優(yōu)勢(shì)更加明顯，故只考慮迭代法的并行化。通過對(duì)迭代法的分析(以穩(wěn)定雙共軛梯度迭代法-DLUSTB為例)，其主要包括三個(gè)功能區(qū)：格式轉(zhuǎn)換、預(yù)處理(LDU分解)和迭代求解等。為了進(jìn)一步了解其在求解模塊中的執(zhí)行時(shí)間，在上述方案中，也將對(duì)應(yīng)的執(zhí)行時(shí)間進(jìn)行了深入分析。格式轉(zhuǎn)換和迭代求解是主要的耗時(shí)部分，而且隨著模型規(guī)模的增大，迭代求解的執(zhí)行時(shí)間會(huì)更長(zhǎng)(圖4)。

圖4 求解模塊中不同功能區(qū)的執(zhí)行時(shí)間Fig.4 Execution time of different functional areas in the solving module

綜上分析，EOS模塊、MULTI模塊、格式轉(zhuǎn)換和迭代求解是算法改進(jìn)和并行化的重點(diǎn)。

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

2.1 OpenMP簡(jiǎn)介

OpenMP是一個(gè)為共享存儲(chǔ)的多處理機(jī)上編寫并行程序而設(shè)計(jì)的編程接口，其功能的實(shí)現(xiàn)得到兩種形式的支持：編譯指導(dǎo)語句和運(yùn)行時(shí)庫函數(shù)，并且可靈活地控制程序運(yùn)行。OpenMP使用的是Fork-Join(交叉-合并)并行執(zhí)行模型。當(dāng)程序運(yùn)行到并行區(qū)域時(shí)，會(huì)創(chuàng)建新線程或者喚醒已有線程(Fork過程)，此時(shí)主線程和派生線程共同工作；當(dāng)執(zhí)行完并行區(qū)域，派生線程會(huì)退出或掛起，回到主線程(Join過程)。OpenMP具有與標(biāo)準(zhǔn)Fortran、C和C++等語言無縫集成，只要在程序中加入OpenMP編譯命令，即可完成并行；通過一些常見函數(shù)設(shè)置、獲取并行環(huán)境參數(shù)，如omp_set_num_threads函數(shù)來設(shè)置程序的線程數(shù)、omp_get_thread_num獲取當(dāng)前線程號(hào)等。

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

基于TOUGH2模擬器算法以及OpenMP并行化的特點(diǎn)，按以下幾方面進(jìn)行了改進(jìn)并實(shí)現(xiàn)了并行。

(1)動(dòng)態(tài)內(nèi)存分配

TOUGH2源程序數(shù)組采用靜態(tài)內(nèi)存分配，導(dǎo)致部分?jǐn)?shù)組擁有大量空間而浪費(fèi)，而部分?jǐn)?shù)組則不夠，帶來十分不合理的內(nèi)存使用方式，同時(shí)對(duì)于大規(guī)模計(jì)算問題，常常無法編譯成功?；诖耍捎妹嫦?qū)ο蟮乃枷牒蛣?dòng)態(tài)內(nèi)存分配的方式重組了源程序。

(2)并行化方案

通過前面對(duì)模擬器中的耗時(shí)分析，程序的計(jì)算時(shí)間主要分布在EOS模塊、MULTI模塊和方程求解等部分。故在進(jìn)行并行化時(shí)，主要針對(duì)這些區(qū)域進(jìn)行并行化，而且對(duì)OpenMP程序本身運(yùn)行的線程鎖、同步開銷和負(fù)載不均衡等關(guān)鍵問題也進(jìn)行了優(yōu)化。

(3)跳轉(zhuǎn)語句

由于OpenMP并行的限制，循環(huán)語句中應(yīng)該是單入口與單出口，不易使用跳轉(zhuǎn)語句，故需對(duì)程序中的跳轉(zhuǎn)語句進(jìn)行處理。其中一部分采用原子鎖來處理，另一部分則需對(duì)程序進(jìn)行重構(gòu)。

(4)參數(shù)優(yōu)化

源程序中各函數(shù)中的變量以及中間變量非常多，若在并行時(shí)都采用私有變量來處理(private、privatefirst、lastprivate或threadprivate)，則會(huì)帶來大量的開銷，對(duì)并行的優(yōu)化效果非常不利。鑒于此，對(duì)程序采取兩種處理方式：第一是將具有多個(gè)功能塊的函數(shù)，分解為多個(gè)獨(dú)立的子函數(shù)，一方面可以減小局部的參數(shù)傳遞，另一方面程序的可讀性也得到了增強(qiáng)。第二是對(duì)多余的參數(shù)進(jìn)行剔除或替換，重要的參數(shù)在函數(shù)之間采用參數(shù)傳遞。

(5)相關(guān)性的處理

一般來說，若兩個(gè)計(jì)算塊之間具有某種相關(guān)性，則它們無法或很難進(jìn)行并行化，如在組建Jacobian矩陣計(jì)算單元之間質(zhì)量與能量項(xiàng)、進(jìn)行的不完全LDU分解、以及稀疏矩陣進(jìn)行乘等部分。這時(shí)可借助臨時(shí)數(shù)組來記錄中間數(shù)據(jù)，后再進(jìn)行統(tǒng)一處理。

(6)并行化原則

在進(jìn)行OpenMP并行化時(shí)，主要遵從三個(gè)原則：第一，減小功能之間或數(shù)據(jù)之間的依賴性，盡可能有更大的并行區(qū)域，減少并行區(qū)域與串行過程的相互交叉；第二，采取更大功能范圍內(nèi)的并行，減小局部功能的并行，以減少并行線程申請(qǐng)的時(shí)間開銷；第三，縮減并行區(qū)域的私有變量個(gè)數(shù)、原子鎖、時(shí)間同步等耗時(shí)操作，并考慮負(fù)載不均衡等問題；第四，盡量保持原來程序的風(fēng)格。

3 并行實(shí)驗(yàn)結(jié)果分析

3.1 模型設(shè)計(jì)及介紹

模型是以鄂爾多斯CCS示范場(chǎng)地的石千峰組地層為研究對(duì)象，模型的水平范圍是以注入井為中心的5 km，厚度為291 m(圖5)。在進(jìn)行并行模擬器測(cè)試時(shí)，主要考慮以下兩個(gè)方面：第一是改進(jìn)后的模擬器與原有模擬器的正確性與效率比較；第二是并行化模擬器的并行化效率。針對(duì)第二方面，設(shè)計(jì)了三個(gè)模型，分別代表小、中與中大規(guī)模，來探討模擬器的并行加速情況，具體模型方案設(shè)計(jì)見表2。

圖5 模型測(cè)試網(wǎng)格Fig.5 Model test grid

序號(hào)模型描述單元連接數(shù)模型1剖面二維模型2 500(50×50)4 900模型2加密模型110 000(100×100)19 800模型3加密模型240 000(200×200)79 600

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

本次實(shí)驗(yàn)在PC機(jī)上完成，CPU為Intel(R) Core(TM) i7-2600 CPU@3.4 GHz，內(nèi)存8 GB，操作系統(tǒng)為32位Windows 7，編譯器采用Microsoft開發(fā)的支持OpenMP2.0的Visual Studio 2005。

3.3 模擬器結(jié)果對(duì)比

在程序編譯時(shí)取消OpenMP支持，則改進(jìn)后的模擬器即為單機(jī)版。下面針對(duì)方案1，分別運(yùn)用改進(jìn)前后的模擬器進(jìn)行運(yùn)算，并進(jìn)行結(jié)果對(duì)比分析。圖6是改進(jìn)前后模擬器運(yùn)算的CO2不同封存形式歷時(shí)曲線。可知改進(jìn)前后模擬器的運(yùn)算結(jié)果穩(wěn)合較好。而且改進(jìn)前后模擬器的運(yùn)算結(jié)果數(shù)值是一樣的(只存在結(jié)果輸出的舍入誤差，表3)。故不難得出，改進(jìn)后的模擬器是可靠的、正確的。

圖6 CO2封存量歷時(shí)曲線對(duì)比(氣、液、總)Fig.6 Duration curves comparison of CO2storage capacity (in gas, liquid, and total)

圖7是改進(jìn)前后兩種模擬器的運(yùn)算總時(shí)間對(duì)比；由圖上的藍(lán)線可知，改進(jìn)后模擬器的運(yùn)算時(shí)間比原模擬器少，最大可達(dá)到1.28倍的加速，最小可達(dá)到1.12倍的加速。

表3 模擬器改進(jìn)前后模擬結(jié)果對(duì)比Table 3 Comparison of simulation results before and after the improvement simulator

圖7 新舊模擬器運(yùn)算時(shí)間對(duì)比(模型1)Fig.7 Comparison of computing time between the new and old simulator in the first model

由以上對(duì)比不難得出，改進(jìn)后模擬器的運(yùn)算結(jié)果是可靠的，而且運(yùn)算效率比原來高。

3.4 并行模擬器模擬結(jié)果對(duì)比

常用的浮點(diǎn)計(jì)算機(jī)都提供符合IEEE-754標(biāo)準(zhǔn)的單精度和雙精度運(yùn)算，常見的計(jì)算機(jī)系統(tǒng)都能夠支持4倍精度運(yùn)算，故在進(jìn)行浮點(diǎn)運(yùn)算時(shí)都會(huì)存在一定的舍入誤差。鑒于此，本文對(duì)不同核數(shù)的模擬器也進(jìn)行了對(duì)比研究，以了解在當(dāng)前并行框架下，并行模擬器的誤差大小和穩(wěn)定性。

圖8是不同核數(shù)對(duì)模型1運(yùn)算的CO2封存量與單機(jī)運(yùn)算CO2封存量之差歷時(shí)曲線；對(duì)比分析有：2核的運(yùn)算結(jié)果與單核的運(yùn)算結(jié)果相差為零；4核與單核的相差最大，最后基本穩(wěn)定在20 kg，這個(gè)誤差相對(duì)于CO2的總封存量(20 a時(shí)，1.84491683E+07 kg)非常小，誤差僅0.000108%。故在當(dāng)前的并行框架下，并行模擬器的運(yùn)算結(jié)果是可靠的、穩(wěn)定的。

圖8 不同核之間CO2封存總量之差Fig.8 The difference between the total amount of CO2 sequestration in different CPU numbers

3.5 加速效率對(duì)比

加速比(Sp)和執(zhí)行效率(E)是測(cè)試并行算法加速效果的2個(gè)重要指標(biāo)。由表4對(duì)比分析可知：并行化處理可明顯加快模擬器的運(yùn)算速度，即加速比會(huì)隨著核數(shù)的增加而增大，原因在于將計(jì)算任務(wù)分解到多個(gè)處理器上，減少了單個(gè)處理器的計(jì)算量，因而可顯著地減小計(jì)算時(shí)間。在當(dāng)前模型中，4核可以達(dá)到2.28倍的加速比。但隨著線程數(shù)的增加，模型的執(zhí)行效率呈現(xiàn)減小的趨勢(shì)；這主要是由于隨著線程數(shù)的增多、線程的創(chuàng)建與銷毀，以及各個(gè)線程之間的負(fù)載不能達(dá)到完全的負(fù)載平衡，由此帶來了時(shí)間延遲而拖慢了執(zhí)行效率。同時(shí)由于模擬器中還存在許多串行部分，并行效率還有待進(jìn)一步提高，以達(dá)到理想的并行效率。

為了進(jìn)一步了解各個(gè)并行塊的執(zhí)行情況，本文也將其執(zhí)行時(shí)間進(jìn)行了對(duì)比分析。圖9是模型2中主要并行塊的執(zhí)行時(shí)間歷時(shí)曲線。EOS模塊和MULTI模塊都可得到明顯的加速；方程求解的加速稍小。通過對(duì)方程求解中不同功能塊進(jìn)行對(duì)比，格式轉(zhuǎn)換的并行加速效果較明顯，預(yù)處理和迭代求解的加速效果稍小。

表4 不同處理核數(shù)下并行處理加速情況Table 4 The parallel processing acceleration in different CPU numbers

圖9 模型2中不同核的執(zhí)行時(shí)間對(duì)比Fig.9 Execution time comparison of different CPU numbers in the second model

4 結(jié)論

(1)通過分析改進(jìn)前后的模型器計(jì)算結(jié)果：改進(jìn)后模擬器的算法是正確可靠的，且改進(jìn)后模擬器的執(zhí)行效率得到了提高，在2 500個(gè)單元的模型中，執(zhí)行速度提高1.12～1.28倍。

(2)由并行模擬器模擬結(jié)果對(duì)比分析，不同核之間由于計(jì)算機(jī)的浮點(diǎn)運(yùn)算誤差而引入的數(shù)值離散是非常小的，故在OpeMP框架下研發(fā)的并行模擬器是可靠的、穩(wěn)定的。分析方案2的運(yùn)算執(zhí)行情況，在4核時(shí)，可以達(dá)到2.28倍的加速。

(3)由于改進(jìn)后的模擬器還存在一些串行執(zhí)行部分，以致未能達(dá)到理想的并行效率；同時(shí)與已有的并行模擬器進(jìn)行對(duì)比研究，將是以后的研究方向。