劉任濤，陳　衛(wèi)

LIU Rentao1，2,CHEN Wei1

1.中國科學(xué)院 計算機(jī)網(wǎng)絡(luò)信息中心 高性能計算技術(shù)與應(yīng)用發(fā)展部，北京 100190

2.中國科學(xué)院大學(xué)，北京 100049

1.Department of High Performance Computing Technology and Application Development,Computer Network Information Center,ChineseAcademy of Sciences,Beijing 100190,China

2.University of ChineseAcademy of Sciences,Beijing 100049,China

1　引言

分治算法是計算機(jī)領(lǐng)域中一種廣為人知的算法，它的功能體現(xiàn)在于兩個方面：其一，通過多次的遞歸迭代自動地將復(fù)雜的問題分解成許許多多小問題；其二，將每個小問題的解沿著分解該問題的逆方向合并，直到得到原問題的解[1]。分治算法實(shí)現(xiàn)這兩個功能的具體方式在于將程序段多次遞歸迭代。近年來，由于能夠以相對簡單的程序段解決復(fù)雜的問題，分治算法的分而治之的思想被運(yùn)用在諸多領(lǐng)域：例如在計算機(jī)科學(xué)領(lǐng)域利用分治算法研究聚類問題和群體智能算法[2-3]；在計算機(jī)圖形學(xué)領(lǐng)域利用分治策略劃分Delaunay三角網(wǎng)格[4-5]；在數(shù)值計算領(lǐng)域利用分治算法求矩陣特征值[6-7]；在系統(tǒng)科學(xué)領(lǐng)域利用分治算法研究種群的進(jìn)化[8]。然而分治策略很少被應(yīng)用于與計算物理和化學(xué)物理相關(guān)的計算機(jī)模擬當(dāng)中。

多孔材料是一種具有非均勻內(nèi)部結(jié)構(gòu)的系統(tǒng)，在物理學(xué)[9-12]、化學(xué)[13]、材料科學(xué)[14]等領(lǐng)域引發(fā)了越來越多研究者的關(guān)注。將多孔材料體系離散為格點(diǎn)系統(tǒng)來進(jìn)行計算機(jī)模擬是一種高效的計算機(jī)模擬方法，在基礎(chǔ)研究領(lǐng)域不斷涌現(xiàn)出新成果[15-16]。多孔材料的內(nèi)部由兩部分組成：一部分是母體，被固體物質(zhì)填充；另一部分是孔隙，可以容納流體（例如氣體和液體）。多孔材料格點(diǎn)模型（Porous Material with Lattice Model，PMLM）的計算機(jī)模擬分為初始化建模和體系演化兩個步驟，其中的初始化建模步驟需要把全體空間中的格點(diǎn)區(qū)分為兩類進(jìn)行標(biāo)記：一類格點(diǎn)標(biāo)記為母體填充物（例如固體）；另一類格點(diǎn)標(biāo)記為材料當(dāng)中的孔隙。然而在實(shí)際研究過程中，對格點(diǎn)進(jìn)行分類的過程并不是完全隨機(jī)的。為了逼近真實(shí)物理體系，通常母體區(qū)塊采用具有硬球勢或Lennard-Jones勢的球體[9，15]。原始建模方法不得不在固定每一個母體球的條件下遍歷地搜索全體空間的每一個格點(diǎn)，其計算量與系統(tǒng)體積的平方成正比。隨著模擬體系規(guī)模的增大，計算機(jī)模擬過程中的初始化建模步驟的計算量急劇增加，這是許多關(guān)于大規(guī)模多孔材料的研究無法承受的。

k近鄰搜索是一種用來取代全空間遍歷搜索的高效的搜索方式，在模式識別、機(jī)器學(xué)習(xí)、數(shù)據(jù)挖掘、時間序列分析等方面有諸多應(yīng)用。利用k近鄰搜索可以預(yù)測數(shù)據(jù)樣本的分類，或以較小的代價尋找目標(biāo)樣本[17-20]。受到k近鄰分類算法的啟發(fā)，本文把基于分治策略的近鄰搜索（Divide-Conquer based Nearest-Neighbor searching，DCNNs）方法應(yīng)用到PMLM的建模當(dāng)中，詳細(xì)介紹了在PMLM系統(tǒng)中應(yīng)用DCNNs建模的方法。本文還在理論和實(shí)驗(yàn)上對原始建模方法和DCNNs建模方法的運(yùn)行效率作了比較。

2　基于分治策略的近鄰搜索方法理論

2.1　PMLM系統(tǒng)的基于分治策略近鄰搜索的優(yōu)化建模方法

在一個三維格點(diǎn)模型中，格點(diǎn)的編號遵循：首先沿x方向布滿再向y方向延拓，布滿每一個xy平面后再向z方向延拓。因此，只要得到一個格點(diǎn)的序號就可以知道該格點(diǎn)在三維體系中的坐標(biāo)(x,y,z)；反之，根據(jù)任意一個格點(diǎn)的坐標(biāo)就可以知道該格點(diǎn)的序號。

多孔材料的顯著特征就是構(gòu)型當(dāng)中的母體物質(zhì)與孔隙分散的占據(jù)全體空間。最常采用的硬球勢體系須將母體區(qū)域視為由一個個小球堆疊而成，孔隙區(qū)域占據(jù)剩余的空間。隨機(jī)選擇一個格點(diǎn)作為一個母體小球的球心A(x0,y0,z0)。為了標(biāo)記出被球形母體占據(jù)的格點(diǎn)，最簡單的篩選辦法就是首先固定母體球的球心坐標(biāo)，然后計算各個格點(diǎn)與球心之間的幾何距離：

球形母體區(qū)域的限制條件：如果某個格點(diǎn)到該球心之間的幾何距離di小于設(shè)定的母體球半徑R，則認(rèn)為此格點(diǎn)在該母體球內(nèi)部。

本文提出的DCNNs建模方法是以PMLM的一個母體區(qū)域的中心A作為根節(jié)點(diǎn)，通過分治的辦法將近鄰區(qū)域的格點(diǎn)劃分為不相交的區(qū)域，遞歸地向外部空間擴(kuò)展以實(shí)現(xiàn)近鄰搜索。利用該DCNNs可以快速地標(biāo)記出PMLM當(dāng)中應(yīng)當(dāng)被母體占據(jù)的格點(diǎn)。

圖1顯示了用于PMLM建模過程中的DCNNs在不同類型格點(diǎn)位置上向近鄰格點(diǎn)擴(kuò)展的狀態(tài)空間圖，圖中的每一個箭頭對應(yīng)于一個方向上的擴(kuò)展（體現(xiàn)在計算機(jī)程序中就是一個子程序）。算法對應(yīng)的偽代碼如算法1～7所示。

圖1　近鄰搜索各個格點(diǎn)的狀態(tài)空間圖

偽代碼算法1表述了向6個方向擴(kuò)展根節(jié)點(diǎn)的算法，狀態(tài)空間如圖1（a）所示。

偽代碼算法2和偽代碼算法3分別表述了向5個方向擴(kuò)展z+節(jié)點(diǎn)或z-節(jié)點(diǎn)的算法，狀態(tài)空間如圖1（b）所示。

算法 2 UP（z，queue）

偽代碼算法4和偽代碼算法5分別表述了向3個方向擴(kuò)展y+節(jié)點(diǎn)或y-節(jié)點(diǎn)的算法，狀態(tài)空間如圖1（c）所示。

DCNNs過程中的數(shù)據(jù)結(jié)構(gòu)可以描繪成樹狀圖，如圖2所示，除了葉節(jié)點(diǎn)以外：

（1）根節(jié)點(diǎn)有6個子節(jié)點(diǎn)；

（2）每個z+或z-節(jié)點(diǎn)有5個子節(jié)點(diǎn)；

（3）每個y+或y-節(jié)點(diǎn)有3個子節(jié)點(diǎn)；

（4）每個x+或x-節(jié)點(diǎn)有1個子節(jié)點(diǎn)。

圖2　三維格點(diǎn)系統(tǒng)的DCNNs的搜索樹

從幾何空間上看，這樣的搜索機(jī)制類似于“碰壁”：設(shè)想在母體球外包裹了一層球殼，當(dāng)計算機(jī)沿著一條搜索路徑向前搜索直到訪問到球殼上的格點(diǎn)時返回尋找下一條路徑繼續(xù)搜索。該DCNNs方式在搜索順序方面類似于深度優(yōu)先搜索，搜索方向沿著搜索樹直到葉節(jié)點(diǎn)才返回，一言以蔽之就是“一路到底”。但區(qū)別在于：傳統(tǒng)的深度優(yōu)先搜索的目的只是為了找到一個目標(biāo)節(jié)點(diǎn)，因此一旦找到目標(biāo)節(jié)點(diǎn)就會終止搜索；而本文介紹的搜索方法并沒有特定的目標(biāo)節(jié)點(diǎn)，只有把搜索樹的所有葉節(jié)點(diǎn)全部搜索了一遍以后才會終止。

該基于分治策略的近鄰搜索方法吸收了數(shù)據(jù)分割的思想[20]，通過將近鄰空間格點(diǎn)劃分成合理的數(shù)據(jù)結(jié)構(gòu)，保證了在同一次近鄰搜索當(dāng)中既不會遺漏任何一個應(yīng)當(dāng)被搜索到的格點(diǎn)，又不會重復(fù)搜索任何一個格點(diǎn)。采用該DCNNs方法對PMLM系統(tǒng)進(jìn)行建模，能夠準(zhǔn)確地篩選出初始構(gòu)型中應(yīng)當(dāng)被母體物質(zhì)占據(jù)的格點(diǎn)。此外，只要改變程序中對于球形母體區(qū)域的限制條件，就可以將該DCNNs方法推廣用來對具有更復(fù)雜的母體形狀的PMLM系統(tǒng)進(jìn)行建模。

2.2　建模時間復(fù)雜度的理論分析與比較

2.2.1原始方法建模的時間復(fù)雜度

對于PMLM系統(tǒng)，原始的算法是：首先隨機(jī)生成每一個母體區(qū)塊的球心的坐標(biāo)，然后在固定該坐標(biāo)的條件下遍歷搜索空間中的所有格點(diǎn)，計算各個格點(diǎn)的坐標(biāo)與母體球心坐標(biāo)的空間距離，從而篩選出位于該母體球內(nèi)部的格點(diǎn)。

不妨設(shè)整個PMLM體系的格點(diǎn)總數(shù)為N；體系中被母體物質(zhì)填充的格點(diǎn)個數(shù)與格點(diǎn)總數(shù)的比例是固定的，設(shè)為η(0＜η＜1)；同一體系中的各個母體球的半徑大小是相同的，設(shè)每個母體球占據(jù)n個格點(diǎn)。對于不同的體系，當(dāng)母體球的半徑大小固定不變時，每個母體球中包含的格點(diǎn)個數(shù)n也是固定的常數(shù)（對于多孔材料n?N）。格點(diǎn)系統(tǒng)建模過程中的時間復(fù)雜度：

2.2.2DCNNs方法建模的時間復(fù)雜度

同樣，不妨設(shè)整個PMLM體系的格點(diǎn)總數(shù)為N；體系中被母體物質(zhì)填充的格點(diǎn)個數(shù)與總的格點(diǎn)個數(shù)的比例為η(0＜η＜1)；每個母體球占據(jù)n個格點(diǎn)（n固定，且n?N）。對于同樣的PMLM系統(tǒng)，基于分治策略的近鄰搜索首先也是隨機(jī)生成每一個母體球的球心的坐標(biāo)，然后在固定該坐標(biāo)的條件下利用2.1節(jié)當(dāng)中介紹的規(guī)則進(jìn)行分類，在每一個母體球球心附近搜索到的格點(diǎn)個數(shù)，僅僅是母體球理論上包含格點(diǎn)個數(shù)n加上母體球外的“壁”包含的格點(diǎn)個數(shù)（其數(shù)量級一般不高于n具有的數(shù)量級）。格點(diǎn)系統(tǒng)建模過程中的時間雜度：

很明顯，在理論上，對于n?N的情形（例如多孔材料），DCNNs方法建模的時間復(fù)雜度更低，能夠節(jié)省大量計算量。

3PMLM系統(tǒng)氣液相變的計算機(jī)模擬

3.1　Gibbs系綜下基于密度泛函的格點(diǎn)模型平均場理論的熱力學(xué)函數(shù)

假定空間被離散化為三維的格點(diǎn)網(wǎng)絡(luò)，共有N個格點(diǎn)，第i個格點(diǎn)的粒子數(shù)密度為ρi，設(shè)第i個格點(diǎn)上的外場為Vi，第i、j格點(diǎn)之間的相互作用強(qiáng)度為εij，則整個體系的總勢能：

對于經(jīng)典體系的熱力學(xué)正則系綜（固定的溫度、體積、總粒子數(shù)），全體粒子所處的密度狀態(tài)服從熱力學(xué)統(tǒng)計規(guī)律，任意一個給定的粒子數(shù)分布{ρi}出現(xiàn)的熱力學(xué)概率為：

其中，正則系綜配分函數(shù)：

因此，體系的Helmhotz自由能：

構(gòu)造一個無外場無相互作用(Vi=0,εij=0)的參考體系，設(shè)參考體系的總能量為E?，則實(shí)際體系的能量E可以寫為參考體系能量與附加能量ΔE的加和：

在參考體系中，各個熱力學(xué)構(gòu)型出現(xiàn)的幾率相等：

此時配分函數(shù)Z的物理意義可以形象地表述為M個粒子隨機(jī)分配給N個格點(diǎn)分配方式總個數(shù)，即服從超幾何分布：

對于M和N的值很大的情形（例如宏觀熱力學(xué)體系），運(yùn)用Stirling近似ln(N!)≈NlnN-N ，則該參考體系的Helmhotz自由能：

其中，kB是Boltzmann常數(shù)，一個關(guān)于溫度和能量的物理常數(shù)。

由于該參考體系的粒子數(shù)分布可以視為宏觀上均勻的，設(shè)參考體系的粒子數(shù)密度ρ=M/N，因而該參考體系的Helmhotz自由能：

其中的密度ρ是無量綱的約化密度(0＜ρ＜1)。

引入平均場近似，設(shè)外場對每一個格點(diǎn)的作用強(qiáng)度為?i，則附加能量可以表示成：

3.2　密度迭代演化的計算機(jī)模擬方法

對于實(shí)際的物理系統(tǒng)，系統(tǒng)達(dá)到平衡時的密度分布使式（11）取得極小值，即滿足微分方程：

將格點(diǎn)之間的相互作用當(dāng)作短程力處理，僅僅考慮最近鄰和次近鄰格點(diǎn)之間的相互作用[14]，并將表達(dá)式（11）代入式（12）即可得到達(dá)到熱力學(xué)平衡狀態(tài)下系統(tǒng)關(guān)于密度的代數(shù)方程[14]：

對于每一個格點(diǎn)，每次把舊的密度值代入上式右端，計算結(jié)果賦值給等號左端的新的密度值ρrhsi。為了防止密度取值跳動過大以致溢出（必須滿足ρ≤1），密度值與舊密度值進(jìn)行線性組合，使得新得到的密度值僅僅在原有值的基礎(chǔ)上小范圍改動。因此，整個迭代的完整過程用等式（14）、（15）、（16）表示[14]：

等式（15）中的系數(shù)α(α?1)用來限制每圈迭代的速率；等式（16）用來標(biāo)定體系的總粒子數(shù)守恒，得到每一輪迭代被更新的密度值。

反復(fù)迭代等式（14）、（15）、（16），使得系統(tǒng)的密度分布最終收斂。當(dāng)某一圈迭代過后系統(tǒng)中密度值的改變量大小的平均值小于設(shè)定的閾值時，認(rèn)為體系已經(jīng)達(dá)到熱力學(xué)近平衡，終止迭代。

3.3　計算機(jī)模擬結(jié)果的比較

本文對采用周期性邊界條件的邊界長為50、65、80、100的立方體形格點(diǎn)模型分別進(jìn)行計算機(jī)模擬。使用原始方法和DCNNs方法分別進(jìn)行建模能夠構(gòu)造相同的初始構(gòu)型（流體的初始密度均勻ρ=0.28），并經(jīng)過迭代演化達(dá)到同樣的熱力學(xué)平衡態(tài)。一個包含106個格點(diǎn)的三維體系的模擬結(jié)果，在x=0的截面的密度分布如圖3所示。圖右側(cè)標(biāo)出了各種顏色代表的流體密度值（無量綱約化密度0＜ρ＜1），圖中空白區(qū)域是建模時建立的母體（流體密度為零），圖中密度稠密的團(tuán)狀聚集區(qū)域是液相凝聚區(qū)，證明該體系已經(jīng)發(fā)生了氣液相變。

圖3　計算機(jī)模擬多孔材料格點(diǎn)模型氣液相變達(dá)到熱力學(xué)近平衡時在x=0截面上的密度分布

盡管使用原始的建模方法與DCNNs建模方法能夠得到相同的物理結(jié)果，但是兩種方法的運(yùn)行時間差異巨大。將程序在初始化和迭代過程中占用的運(yùn)行時間進(jìn)行比較，如圖4所示。隨著格點(diǎn)系統(tǒng)體積的增加，用DCNNs方法進(jìn)行建模初始化占用的時間是線性增加的，系統(tǒng)遵循熱力學(xué)規(guī)律的迭代計算時間也近似于線性遞增，然而用原始的建模方法進(jìn)行初始化占用的運(yùn)行時間隨著系統(tǒng)的格點(diǎn)總數(shù)的二次方急劇增加。

圖4　不同規(guī)模體系模擬的各步驟的運(yùn)行時間

當(dāng)系統(tǒng)體積達(dá)到106個格點(diǎn)時，應(yīng)用原始的建模方法與DCNNs建模方法進(jìn)行計算機(jī)模擬的過程中各步驟占用時間的比例，如圖5所示。

圖5　N=106的體系在計算機(jī)模擬過程中運(yùn)行時間的組成

不難發(fā)現(xiàn)，隨著體系規(guī)模的增大，若采用原始建模方法則初始化步驟會占據(jù)大部分的運(yùn)行時間。加之系統(tǒng)演化過程可以很方便地并行加速，初始化步驟的運(yùn)行時間在總運(yùn)行時間當(dāng)中占據(jù)的比例有可能會更大。因而對于采用原始方法建模的大規(guī)模PMLM系統(tǒng)，初始化步驟嚴(yán)重拖累了整個計算機(jī)模擬過程的運(yùn)行時間，成為用計算機(jī)模擬大規(guī)模PMLM系統(tǒng)的技術(shù)瓶頸。應(yīng)用DCNNs方法進(jìn)行建模則只需占用幾乎可以忽略不計的運(yùn)行時間來完成初始化，為更大規(guī)模PMLM系統(tǒng)的計算機(jī)模擬創(chuàng)造了條件。

4　結(jié)束語

本文主要研究了一種適用于PMLM系統(tǒng)進(jìn)行計算機(jī)建模的優(yōu)化方法。將分治思想與近鄰搜索策略相結(jié)合，提出的基于分治策略的近鄰搜索方法能將計算機(jī)建模PMLM系統(tǒng)的時間復(fù)雜度由O(N2)降低至O(N)。對PMLM系統(tǒng)氣液相變過程的計算機(jī)模擬，驗(yàn)證了該方法在不影響物理結(jié)果的前提下的高效性。本文提出的優(yōu)化方法突破了計算機(jī)模擬大規(guī)模PMLM系統(tǒng)的一個技術(shù)瓶頸。

參考文獻(xiàn)：

[1]Cormen T H，Leiserson C E，Rivest R L，et al.算法導(dǎo)論：Introduction to algorithms[M].殷建平，徐云，王剛，等譯.3版.北京：機(jī)械工業(yè)出版社，2013：16-17.

[2]王寶文，閻俊梅，劉文遠(yuǎn)，等.基于分治法的高維大數(shù)據(jù)集模糊聚類算法[J].計算機(jī)工程，2007，33（24）：60-62.

7月中旬，我們形成了10個文件，經(jīng)浦東開發(fā)領(lǐng)導(dǎo)小組、市政府常務(wù)會、市委常委會和市人大常委會相繼審議通過。

[3]李田來，劉方愛，王新華.基于分治策略的改進(jìn)人工蜂群算法[J].控制與決策，2015，30（2）：316-320.

[4]余杰，呂品，鄭昌文.Delaunay三角網(wǎng)構(gòu)建方法比較研究[J].中國圖象圖形學(xué)報，2010，15（8）：1158-1167.

[5]Cignonit P，Montanit C，Scopigno R.DeWall：A fast divide and conquer delaunay triangulation algorithm in Ed[J].Computer-Aided Design，1998，30（5）：333-341.

[6]Gansterer W N，Ward R C，Muller R P.An extension of the divide-and-conquer method for a cla-ss of symmetric block-tridiagonaleigenproblems[J].ACMTransactions on Mathematical Software，2002，28（1）：45-58.

[7]魏立峰，李曉梅.求解對稱帶狀廣義特征值問題的擴(kuò)展分治算法[J].計算機(jī)研究與發(fā)展，2004，41（5）：861-867.

[8]王攀，萬君康，馮珊，等.一類基于分治原理的多種群協(xié)同進(jìn)化算法[J].系統(tǒng)工程與電子技術(shù)，2004，26（11）：1687-1690.

[9]Brennan J K，Dong W.Phase transition of one-component fluids absorbed in random porous media：Monte Carlo simulations[J].The Journal of Chemical Physics，2002，116（20）：8948-8958.

[11]馬強(qiáng)，陳振乾.分形多孔材料雙尺度孔隙內(nèi)氣體脫附擴(kuò)散過程數(shù)值模擬[J].東南大學(xué)學(xué)報：自然科學(xué)版，2015，45（4）：728-733.

[12]劉高潔，郭照立，施保昌.多孔介質(zhì)中流體流動及擴(kuò)散的耦合格子 Boltzmann 模型[J].物理學(xué)報，2016，65（1）：282-290.

[13]馬強(qiáng)，陳俊，陳振乾.分形多孔介質(zhì)傳熱傳質(zhì)過程的格子Boltzmann模擬[J].化工學(xué)報，2014，65（S1）：180-187.

[14]崔靜潔，何文，廖世軍，等.多孔材料的孔結(jié)構(gòu)表征及其分析[J].材料導(dǎo)報，2009，23（7）：82-86.

[15]Hughes A P，Thiele U，Archer A J.An introduction to inhomogeneousliquids，densityfunctionaltheory，and the wetting transition[J].American Journal of Physics，2014，82：1119-1129.

[16]Dickey J M，Paskin A.Computer simulation of the lattice dynamics of solids[J].Physical Review，1969，188（3）：1407-1418.

[17]祝繼華，尹俊，邗汶鋅，等.面向低維點(diǎn)集配準(zhǔn)的高效最近鄰搜索法[J].模式識別與人工智能，2014，27（12）：1071-1077.

[18]譚駿祥，李少達(dá)，楊容浩.迭代最近點(diǎn)匹配算法的樹結(jié)構(gòu)k近鄰搜索比較研究[J].測繪科學(xué)，2014，39（4）：152-155.

[19]衛(wèi)煒，張麗艷，周來水.一種快速搜索海量數(shù)據(jù)集K-近鄰空間球算法[J].航空學(xué)報，2006，27（5）：944-948.

[20]肇瑩，劉紅星，王仲宇，等.最近鄰搜索用于分類問題的一種改進(jìn)[J].南京大學(xué)學(xué)報：自然科學(xué)版，2009，45（4）：455-462.

[21]Beckmann N，Krigel H P，Schneider R，et al.The R*-tree：An efficient and robust access method for points and rectangles[C]//Hector G M，Jagadish H V.Proceedings of the 1990 ACM SIGMOD Conference on Management of Data.New York：ACM，1990：322-331.