劉任濤，陳　衛(wèi)

LIU Rentao1，2,CHEN Wei1

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

2.中國科學院大學，北京 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　引言

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

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

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

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

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

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

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

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

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

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

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

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

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

算法 2 UP（z，queue）

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

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

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

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

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

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

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

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

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

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

2.2.1原始方法建模的時間復雜度

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

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

2.2.2DCNNs方法建模的時間復雜度

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

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

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

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

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

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

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

因此，體系的Helmhotz自由能：

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

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

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

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

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

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

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

引入平均場近似，設外場對每一個格點的作用強度為?i，則附加能量可以表示成：

3.2　密度迭代演化的計算機模擬方法

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

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

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

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

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

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

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

圖3　計算機模擬多孔材料格點模型氣液相變達到熱力學近平衡時在x=0截面上的密度分布

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

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

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

圖5　N=106的體系在計算機模擬過程中運行時間的組成

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

4　結(jié)束語

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

參考文獻：

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

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

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

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

[4]余杰，呂品，鄭昌文.Delaunay三角網(wǎng)構(gòu)建方法比較研究[J].中國圖象圖形學報，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]魏立峰，李曉梅.求解對稱帶狀廣義特征值問題的擴展分治算法[J].計算機研究與發(fā)展，2004，41（5）：861-867.

[8]王攀，萬君康，馮珊，等.一類基于分治原理的多種群協(xié)同進化算法[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]馬強，陳振乾.分形多孔材料雙尺度孔隙內(nèi)氣體脫附擴散過程數(shù)值模擬[J].東南大學學報：自然科學版，2015，45（4）：728-733.

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

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

[14]崔靜潔，何文，廖世軍，等.多孔材料的孔結(jié)構(gòu)表征及其分析[J].材料導報，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]祝繼華，尹俊，邗汶鋅，等.面向低維點集配準的高效最近鄰搜索法[J].模式識別與人工智能，2014，27（12）：1071-1077.

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

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

[20]肇瑩，劉紅星，王仲宇，等.最近鄰搜索用于分類問題的一種改進[J].南京大學學報：自然科學版，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.