王 姝，王小鴿，楊廣文

清華大學(xué) 計(jì)算機(jī)系，北京100084

1 引言

隨著燃燒數(shù)值模擬計(jì)算技術(shù)的發(fā)展，模擬的問題越來越復(fù)雜，模擬的精度要求越來越高，傳統(tǒng)模擬方法采用結(jié)構(gòu)網(wǎng)格，計(jì)算較為簡單，但是受到其結(jié)構(gòu)特性的限制，不能靈活地根據(jù)模擬情況改變網(wǎng)格的結(jié)構(gòu)以適應(yīng)新的計(jì)算需求，因此逐漸被非結(jié)構(gòu)網(wǎng)格所替代。相對于結(jié)構(gòu)網(wǎng)格而言，非結(jié)構(gòu)網(wǎng)格可以對復(fù)雜的幾何外形結(jié)構(gòu)進(jìn)行更為準(zhǔn)確的描述；非結(jié)構(gòu)網(wǎng)格中的一個(gè)點(diǎn)周圍的點(diǎn)數(shù)和單元數(shù)不是固定的，這一特點(diǎn)使得非結(jié)構(gòu)網(wǎng)格易于調(diào)整網(wǎng)格結(jié)構(gòu)，以適應(yīng)計(jì)算的需要。

在進(jìn)行燃燒數(shù)值模擬計(jì)算的過程中，部分區(qū)域的計(jì)算精度要求比其他區(qū)域高，如化學(xué)反應(yīng)劇烈的區(qū)域、流場梯度變化劇烈的區(qū)域。為了提高計(jì)算精度，需要對此類區(qū)域的網(wǎng)格進(jìn)行自適應(yīng)加密，通過減小計(jì)算網(wǎng)格的尺度以達(dá)到需要的精度。另外，在某些不需要高計(jì)算精度的區(qū)域，為了減少計(jì)算量，對網(wǎng)格進(jìn)行自適應(yīng)減疏。

由于數(shù)值計(jì)算規(guī)模越來越大，傳統(tǒng)的串行程序已經(jīng)不能滿足高效的計(jì)算需求，將數(shù)值模擬計(jì)算并行化成為迫切需求。同時(shí)，在計(jì)算并行化以后，由于計(jì)算集中在化學(xué)反應(yīng)劇烈的區(qū)域，進(jìn)程之間的計(jì)算負(fù)載會嚴(yán)重不平衡，將導(dǎo)致并行效率大大降低。非結(jié)構(gòu)網(wǎng)格和網(wǎng)格自適應(yīng)結(jié)合使用可以達(dá)到模擬計(jì)算的精度要求，但是含自適應(yīng)功能以后，網(wǎng)格數(shù)目的變動(dòng)會使負(fù)載隨之變化。在這種情況下，需要在進(jìn)程之間進(jìn)行負(fù)載平衡調(diào)度，使計(jì)算負(fù)載均勻分布在所有進(jìn)程上，使得非結(jié)構(gòu)網(wǎng)格自適應(yīng)計(jì)算模擬在提高精度的同時(shí)，縮短程序運(yùn)行的時(shí)間，提高數(shù)值模擬計(jì)算的效率。

在自適應(yīng)非結(jié)構(gòu)網(wǎng)格上實(shí)現(xiàn)程序并行及負(fù)載平衡，最大的挑戰(zhàn)是如何在原有數(shù)據(jù)結(jié)構(gòu)的基礎(chǔ)上正確、高效地執(zhí)行程序。本文是基于文獻(xiàn)[1]中的燃燒數(shù)值模擬和文獻(xiàn)[2]中的非結(jié)構(gòu)網(wǎng)格自適應(yīng)方式，在充分研究了非結(jié)構(gòu)網(wǎng)格的數(shù)據(jù)結(jié)構(gòu)和非結(jié)構(gòu)網(wǎng)格自適應(yīng)變化特點(diǎn)的基礎(chǔ)上，提出了一套有效的負(fù)載平衡策略，使得進(jìn)程之間的負(fù)載較為平均，減少了程序總體運(yùn)行時(shí)間；為了使負(fù)載轉(zhuǎn)移不受非結(jié)構(gòu)網(wǎng)格存儲方式的限制，本文對非結(jié)構(gòu)網(wǎng)格的存儲方式進(jìn)行了改進(jìn)，使得負(fù)載轉(zhuǎn)移更加容易實(shí)現(xiàn)。

2 相關(guān)工作

在非結(jié)構(gòu)網(wǎng)格上實(shí)現(xiàn)并行和負(fù)載平衡計(jì)算逐漸成為數(shù)值模擬計(jì)算的一個(gè)重要研究方向，近幾十年有許多這方面的系統(tǒng)和算法。國際上的許多靜態(tài)負(fù)載平衡方法[3-6]是在編譯階段預(yù)測負(fù)載，并將數(shù)據(jù)和計(jì)算映射到進(jìn)程上，這對于常規(guī)并行問題是很有效的，但是當(dāng)問題變得非常規(guī)時(shí)（例如需要引用非本地?cái)?shù)據(jù)或者計(jì)算迭代的開銷在運(yùn)行時(shí)無法預(yù)測）性能會變得很差。SUPPLE[7]結(jié)合了靜態(tài)和動(dòng)態(tài)負(fù)載平衡的方法，首先通過靜態(tài)的方式在編譯階段預(yù)測負(fù)載和通信，程序運(yùn)行時(shí)采用塊為單位進(jìn)行負(fù)載轉(zhuǎn)移。還有很多對負(fù)載的劃分基于圖劃分的方法[8-11]，它們的網(wǎng)格劃分方法受到非結(jié)構(gòu)網(wǎng)格和自適應(yīng)機(jī)制的約束，很難應(yīng)用到現(xiàn)有的數(shù)值模擬計(jì)算中。

在國內(nèi)的研究中，DSMC[12]并行算法的負(fù)載平衡的指標(biāo)是每個(gè)進(jìn)程上的網(wǎng)格上的分子數(shù)量大致相當(dāng)，采用交替二分法劃分網(wǎng)格，這種劃分方法的負(fù)載平衡效率較好，但會產(chǎn)生許多不規(guī)則的通信，容易造成程序效率低下。文獻(xiàn)[13]采用預(yù)處理方法采用METIS 進(jìn)行負(fù)載平衡和任務(wù)劃分，并提出一種網(wǎng)格單元重排序算法提高訪存命中率難點(diǎn)以及研究的目標(biāo)。這種動(dòng)態(tài)負(fù)載平衡效果較好，但是需要改變原有數(shù)據(jù)結(jié)構(gòu)的存儲方式。還有文獻(xiàn)[14]采用子區(qū)域數(shù)遠(yuǎn)遠(yuǎn)小于處理器數(shù)目的方式實(shí)現(xiàn)動(dòng)態(tài)負(fù)載，這樣的處理方式可以得到較為均衡的負(fù)載，但是由于問題劃分得過小，需要經(jīng)常通信，會降低程序執(zhí)行的效率。文獻(xiàn)[15]在全局范圍內(nèi)實(shí)現(xiàn)網(wǎng)格的重新劃分，這種處理方式會帶來大量、經(jīng)常的通信。

國外的負(fù)載平衡系統(tǒng)和算法絕大部分是不開源的，無法應(yīng)用到現(xiàn)有的數(shù)值模擬計(jì)算中。而國內(nèi)的處理方法大都需要靈活網(wǎng)格的存儲方式，例如文獻(xiàn)[13]的網(wǎng)格重排序和劃分、文獻(xiàn)[15]的網(wǎng)格重新劃分等。要在盡量不改變非結(jié)構(gòu)網(wǎng)格自適應(yīng)存儲方式的前提下達(dá)到負(fù)載平衡的目的，需要一種對數(shù)據(jù)存儲方式要求不高的負(fù)載平衡方法。

在自適應(yīng)非結(jié)構(gòu)網(wǎng)格上實(shí)現(xiàn)動(dòng)態(tài)負(fù)載平衡需要考慮其數(shù)據(jù)結(jié)構(gòu)和自適應(yīng)方式。每個(gè)網(wǎng)格單元上的數(shù)值需要由其相鄰網(wǎng)格單元上的數(shù)值計(jì)算得出，如圖1 的偽代碼所示。

圖1 單元格計(jì)算的偽代碼

非結(jié)構(gòu)網(wǎng)格中，相鄰網(wǎng)格的存儲地址是通過兩個(gè)網(wǎng)格共享邊的單元格屬性數(shù)組獲得的。如圖2 所示。

圖2 邊的編號、方向和單元格的關(guān)系

從圖2 中可以看出，編號為I的單元格的四條邊從右邊開始逆時(shí)針編號，依次為NEG1(I) ，NEG2(I)，NEG3(I)，NEG4(I)。邊被賦予了方向，其中水平邊以向右為正方向，豎直邊以向上為正方向。沿著邊的正方向，左邊的單元格為該邊的左單元格，右邊的單元格為該邊的右單元格。編號為J的邊的左右單元格分別為NL0(J)和NR0(J)。要想獲得單元格I的鄰居，以右鄰居為例，則應(yīng)先取它的右邊NEG1(I)，然后取這條右邊的右屬性單元NR0(NEG1(I))，從而得到它的右鄰居。

網(wǎng)格的自適應(yīng)加密方式如圖3 所示，1 號網(wǎng)格加密變成4 個(gè)子單元格，子單元格還可以再加密變成更小的4 個(gè)子單元格，以3 號網(wǎng)格為例。

圖3 網(wǎng)格自適應(yīng)加密方式

偏微分方程求解和網(wǎng)格自適應(yīng)需要用到鄰居單元格的數(shù)據(jù)，而鄰居單元格通過圖2 所示的關(guān)系得到。并行求解時(shí)，兩個(gè)相鄰網(wǎng)格被分配到不同的進(jìn)程上，它們共享的邊（即重疊邊）需要分別在這兩個(gè)進(jìn)程上存放。在自適應(yīng)加密或減疏以后，由于兩個(gè)重疊邊的分裂情況不同，需要更新分裂后子邊對應(yīng)的新重疊邊。如圖4（a）所示，初始時(shí)有兩條均未分裂的重疊邊I，J，邊J先分裂成兩個(gè)子邊，而I未分裂（如圖4（b）所示），這時(shí)邊I的重疊邊仍為J，邊J的子邊SUB1(J)和SUB2(J)的重疊邊為邊I。當(dāng)邊I也發(fā)生分裂時(shí)，邊J子邊的重疊邊才發(fā)生更新。如圖4（c）所示，I的子邊SUB1(I)的重疊邊為SUB1(J)，SUB2(I)的重疊邊為SUB2(J)。

本文使用的實(shí)驗(yàn)是一個(gè)化學(xué)反應(yīng)流計(jì)算程序。程序執(zhí)行1 280 步，采用8 個(gè)進(jìn)程（便于在圖上顯示）。初始網(wǎng)格為16×4，化學(xué)反應(yīng)穩(wěn)定后的網(wǎng)格規(guī)模在20 000 左右?；瘜W(xué)反應(yīng)使用常微分方程求解，流動(dòng)使用偏微分方程求解。通過實(shí)驗(yàn)得出，各進(jìn)程解常微分方程的時(shí)間差異很大，而解偏微分方程的時(shí)間差異不大，因此計(jì)算不平衡部分集中在解常微分方程部分。本文關(guān)注的重點(diǎn)在于求解常微分方程計(jì)算時(shí)間的差異，而非通常意義下的負(fù)載，即在求解過程中做了多少、什么種類的運(yùn)算。為說明簡便，本文后續(xù)內(nèi)容中的負(fù)載均指求解常微分方程的計(jì)算負(fù)載?？紤]進(jìn)程間的負(fù)載變動(dòng)情況，統(tǒng)計(jì)程序執(zhí)行計(jì)算的時(shí)間（s），為使結(jié)果更清晰，將時(shí)間t取對數(shù)ln t，各進(jìn)程的執(zhí)行時(shí)間情況如圖5所示。

圖4 原網(wǎng)格相鄰網(wǎng)格被分配到不同的進(jìn)程示意圖

圖5 不同時(shí)間步上各進(jìn)程的負(fù)載分布

從圖5 可以看出，在程序執(zhí)行初期，從1 號進(jìn)程開始，負(fù)載逐漸變大，并且隨著時(shí)間的推進(jìn)各進(jìn)程的負(fù)載陸續(xù)增加，這是第一階段；直到第400 步左右，7，8 號進(jìn)程的負(fù)載劇烈增加，這是第二階段；到第600 步以后，各進(jìn)程的負(fù)載趨于穩(wěn)定，這是第三階段。出現(xiàn)這種情況的原因是：第一階段化學(xué)反應(yīng)不明顯，流場梯度變化劇烈的區(qū)域網(wǎng)格發(fā)生加密，隨著該區(qū)域的位置不斷轉(zhuǎn)移，各進(jìn)程的負(fù)載隨之變化。在第二階段，在7，8 進(jìn)程的計(jì)算區(qū)域突然發(fā)生了劇烈的化學(xué)反應(yīng)，導(dǎo)致這些區(qū)域的負(fù)載急劇增加。在第三階段，化學(xué)反應(yīng)趨于穩(wěn)定，負(fù)載變動(dòng)緩慢。

在文獻(xiàn)[1]的工作基礎(chǔ)上，為了提高計(jì)算效率，將原有并行程序進(jìn)行優(yōu)化。本文在原有數(shù)據(jù)結(jié)構(gòu)上，根據(jù)非結(jié)構(gòu)網(wǎng)格和網(wǎng)格自適應(yīng)方式的特點(diǎn)，對非結(jié)構(gòu)網(wǎng)格的存儲方式進(jìn)行了優(yōu)化，提出了一套高效的負(fù)載轉(zhuǎn)移方法。

3 需要考慮的問題

要在自適應(yīng)非結(jié)構(gòu)網(wǎng)格上實(shí)現(xiàn)動(dòng)態(tài)負(fù)載平衡，需要考慮以下六個(gè)方面：

（1）非結(jié)構(gòu)網(wǎng)格

在傳統(tǒng)的結(jié)構(gòu)網(wǎng)格上，網(wǎng)格單元使用多維數(shù)組存放，并且可以根據(jù)該數(shù)組的下標(biāo)關(guān)系來訪問鄰接單元格，每個(gè)單元格的鄰接單元格是確定的，這便于網(wǎng)格的剖分；而非結(jié)構(gòu)網(wǎng)格通過一維數(shù)組存放單元格，鄰居單元格的數(shù)據(jù)調(diào)用方式如圖2 所示。

（2）自適應(yīng)

當(dāng)一個(gè)網(wǎng)格單元與其相鄰網(wǎng)格的梯度超過閾值時(shí)，網(wǎng)格在型心點(diǎn)處分裂成四個(gè)，稱為加密。原來的網(wǎng)格稱為父網(wǎng)格，四個(gè)新產(chǎn)生的網(wǎng)格稱為子網(wǎng)格。當(dāng)一個(gè)父網(wǎng)格單元與其相鄰網(wǎng)格的梯度低于某一閾值時(shí)，其四個(gè)子網(wǎng)格合并為一個(gè)，稱為減疏。網(wǎng)格的自適應(yīng)會導(dǎo)致某些網(wǎng)格上的計(jì)算急劇增加，加劇了進(jìn)程之間的負(fù)載不平衡性，更加大了負(fù)載平衡的難度。非結(jié)構(gòu)網(wǎng)格的這種特點(diǎn)使得動(dòng)態(tài)負(fù)載過程中的負(fù)載劃分變得困難。

（3）負(fù)載的變動(dòng)

數(shù)值模擬計(jì)算一般包括兩個(gè)方面：求解常微分方程和求解偏微分方程，本文使用的程序中，化學(xué)反應(yīng)使用常微分方程求解，流動(dòng)使用偏微分方程求解。通過分析程序可以得出：不同進(jìn)程之間的化學(xué)反應(yīng)計(jì)算負(fù)載差異巨大；流動(dòng)計(jì)算僅與網(wǎng)格數(shù)目有關(guān)，而化學(xué)反應(yīng)穩(wěn)定以后，各個(gè)進(jìn)程中的網(wǎng)格數(shù)量基本一致。常微分方程通過牛頓迭代求解，在劇烈化學(xué)反應(yīng)開始之前，常微分方程求解計(jì)算負(fù)載僅與網(wǎng)格數(shù)有關(guān)；在劇烈化學(xué)反應(yīng)開始之后，牛頓迭代次數(shù)急劇上升，同一單元格上的計(jì)算負(fù)載相差幾百倍，因此可以得出結(jié)論：同一個(gè)網(wǎng)格的計(jì)算負(fù)載與當(dāng)前牛頓迭代次數(shù)相關(guān)。

由于負(fù)載的變動(dòng)，需要解決的問題包括負(fù)載預(yù)測的方法和負(fù)載轉(zhuǎn)移的時(shí)機(jī)。即負(fù)載的變動(dòng)與哪些因素有關(guān)，如何將這些因素與負(fù)載聯(lián)系起來，以及在執(zhí)行哪些計(jì)算時(shí)，需要負(fù)載轉(zhuǎn)移。

（4）負(fù)載轉(zhuǎn)移的粒度

根據(jù)數(shù)值計(jì)算系統(tǒng)的計(jì)算特點(diǎn)，求解偏微分方程需要由其鄰居單元格上的數(shù)值計(jì)算得出，采用原始粗網(wǎng)格中的一列作為負(fù)載轉(zhuǎn)移的單位，可以降低負(fù)載轉(zhuǎn)移后尋找鄰居單元格的難度。

（5）非結(jié)構(gòu)網(wǎng)格自適應(yīng)存儲方式帶來的困難

圖6 原始網(wǎng)格

通過結(jié)合圖6 和圖7 來分析非結(jié)構(gòu)網(wǎng)格自適應(yīng)時(shí)存儲數(shù)組的特點(diǎn)。假設(shè)原始網(wǎng)格中有1～9 號非父（沒有子單元格）的單元格，按列先序一起存儲在圖7（a）所示的數(shù)組中。非父單元格占據(jù)數(shù)組的前半部分位置，父單元格占據(jù)數(shù)組的后半部分位置，中間未填滿的部分留空。原始未加密網(wǎng)格的存儲方式，如圖7（a）所示，此時(shí)通過計(jì)算，2 號單元格需要加密。那么2 號單元格的子單元格編號為10，11，12，13，放入非父單元格末尾，如圖7（b）所示。此時(shí)2 號單元格已經(jīng)加密成為父單元格，需要將它挪入父單元格存儲區(qū)域，即從數(shù)組末尾開始往前存放。然后將末尾的13 號單元格填補(bǔ)上2號單元格在非父單元格存儲區(qū)域中留下的空白，如圖7（c）所示。

圖7 網(wǎng)格自適應(yīng)存儲方式

由于負(fù)載轉(zhuǎn)移的單位為以原始粗網(wǎng)格中的一列，而非結(jié)構(gòu)網(wǎng)格的存儲方式是不以列為單位的，這樣一來，要從如圖7 所示的非結(jié)構(gòu)網(wǎng)格存放的一維數(shù)組中挑選出第一列（1，3，10，11，12，13 號單元格）上的數(shù)據(jù)執(zhí)行轉(zhuǎn)移，不便于處理。這里需要解決的問題是如何改進(jìn)網(wǎng)格自適應(yīng)存儲方式，達(dá)到第4 章中提出的以列為單位執(zhí)行負(fù)載轉(zhuǎn)移。

（6）負(fù)載轉(zhuǎn)移策略

負(fù)載轉(zhuǎn)移策略的目標(biāo)是使通信量盡可能少，每個(gè)進(jìn)程上的負(fù)載趨于一致。因此在每次負(fù)載轉(zhuǎn)移之前，需要一個(gè)進(jìn)程統(tǒng)計(jì)所有進(jìn)程的負(fù)載情況，各進(jìn)程負(fù)載的最大值與平均負(fù)載的比值為Max/Avg 。在負(fù)載平衡的情況下，Max/Avg為1。不平衡率越大，負(fù)載越不平衡。在進(jìn)行負(fù)載調(diào)度時(shí)，判斷全局的Max/Avg 是否大于某個(gè)給定的閾值，如果大于，則執(zhí)行負(fù)載轉(zhuǎn)移。這個(gè)閾值的大小是根據(jù)負(fù)載轉(zhuǎn)移的開銷和不平衡情況綜合考慮的。遷移的開銷越大，閾值的設(shè)置也需要大些，以避免負(fù)載平衡工作的“得不償失”。

4 解決方案

本文的實(shí)驗(yàn)使用的計(jì)算節(jié)點(diǎn)采用兩個(gè)Intel Xeon X5670六核處理器（2.93 GHz，12 MB Cache），160 GB SATA 硬盤。節(jié)點(diǎn)內(nèi)存48 GB，結(jié)果取多次實(shí)驗(yàn)的平均值。

結(jié)合上一章中提出非結(jié)構(gòu)網(wǎng)格的特點(diǎn)和網(wǎng)格自適應(yīng)方式，針對需要解決的問題，本文提出的負(fù)載平衡方案具體如下。

4.1 負(fù)載預(yù)測方法

經(jīng)過分析串行程序各部分運(yùn)行時(shí)間，解常微分方程時(shí)間占解常微分方程和解偏微分方程總時(shí)間的34.4%。偏微分方程計(jì)算時(shí)間只與網(wǎng)格數(shù)目有關(guān)，常微分方程計(jì)算約90%的時(shí)間花費(fèi)在牛頓迭代上，且臨近時(shí)間步的牛頓迭代次數(shù)相近，因此將本時(shí)間步牛頓迭代次數(shù)作為下一時(shí)間步負(fù)載的估計(jì)。

為了更加精確地預(yù)測負(fù)載，需要選取能夠進(jìn)一步表征常微分方程計(jì)算負(fù)載的變量進(jìn)行預(yù)測。在本實(shí)驗(yàn)情況下，選取化學(xué)反應(yīng)放熱量Q。在放熱量達(dá)到某個(gè)閾值時(shí)，開始劇烈化學(xué)反應(yīng)，將此閾值作為切換負(fù)載預(yù)測策略的分割點(diǎn)。采用這種自動(dòng)切換負(fù)載預(yù)測方法的原因是：開始劇烈化學(xué)反應(yīng)之前，放熱量Q 不穩(wěn)定，無法作為負(fù)載預(yù)測的判據(jù)，此時(shí)解常微分方程的負(fù)載與牛頓迭代次數(shù)成正比；在劇烈化學(xué)反應(yīng)開始之后，放熱量Q 與負(fù)載的相關(guān)性更高。通過結(jié)合兩種預(yù)測方式，可以使負(fù)載的預(yù)測更接近真實(shí)情況。這種負(fù)載預(yù)測方法的效果在第5 章中分析。

4.2 負(fù)載轉(zhuǎn)移時(shí)機(jī)

考慮到第3 章中提出的負(fù)載變動(dòng)方式，由于解偏微分方程的計(jì)算負(fù)載只與網(wǎng)格數(shù)有關(guān)，化學(xué)反應(yīng)穩(wěn)定后各進(jìn)程上的網(wǎng)格數(shù)目差異不大，求解偏微分方程計(jì)算負(fù)載基本平衡，因此只需在解常微分方程時(shí)執(zhí)行動(dòng)態(tài)負(fù)載平衡策略。

由主控進(jìn)程統(tǒng)計(jì)所有計(jì)算進(jìn)程的負(fù)載，若Max/Avg 大于某個(gè)閾值，則需要進(jìn)行負(fù)載轉(zhuǎn)移。

另外，為了避免自適應(yīng)時(shí)重疊邊的更新成為難點(diǎn)，采取在自適應(yīng)之后執(zhí)行負(fù)載轉(zhuǎn)移這種方式。這樣僅需把計(jì)算負(fù)載轉(zhuǎn)移出去，計(jì)算完成之后更新單元格的狀態(tài)即可。

4.3 數(shù)據(jù)存儲方式的調(diào)整

由第3 章的描述可以得知，非結(jié)構(gòu)網(wǎng)格的存儲結(jié)構(gòu)特點(diǎn)是：單元格在存儲時(shí)沒有規(guī)律可循，要從所有網(wǎng)格中挑出某一列單元轉(zhuǎn)移到其他進(jìn)程上，需要使用數(shù)組標(biāo)記每個(gè)單元格所在的列，當(dāng)執(zhí)行負(fù)載轉(zhuǎn)移時(shí)，遍歷標(biāo)記數(shù)組，從中挑出所需的單元格，并放入一塊連續(xù)的數(shù)組中。為了避免這種方式帶來的麻煩，改為修改非結(jié)構(gòu)網(wǎng)格的存儲方式，將整個(gè)網(wǎng)格存儲地址以列為單位進(jìn)行劃分，如圖8 所示，各列分別占據(jù)一塊連續(xù)的空間，與圖7 相同，2 號單元格加密，其所在列的任何單元格的增加和變動(dòng)都在這一塊空間中，要轉(zhuǎn)移第一列（1，3，10，11，12，13 號單元格）不需挑選，直接轉(zhuǎn)移即可。這使得以列為單位的負(fù)載轉(zhuǎn)移成為可能，并且程序在執(zhí)行時(shí)有較高的空間局部性。

圖8 以列為單位的存儲方式

4.4 負(fù)載轉(zhuǎn)移策略

本文采用的轉(zhuǎn)移策略為進(jìn)程一對一進(jìn)行負(fù)載轉(zhuǎn)移，這種遷移是一種簡單易行的方法，在其基礎(chǔ)上可以實(shí)現(xiàn)更為復(fù)雜的遷移策略。

負(fù)載平衡分為兩個(gè)階段：

負(fù)載轉(zhuǎn)移過程：將進(jìn)程按照預(yù)測的負(fù)載大小排序，負(fù)載轉(zhuǎn)移的規(guī)則為最大負(fù)載的進(jìn)程轉(zhuǎn)移給最小負(fù)載的進(jìn)程，次大負(fù)載的進(jìn)程轉(zhuǎn)移給次小負(fù)載的進(jìn)程，以此類推。在過載進(jìn)程（即進(jìn)程負(fù)載大于平均負(fù)載）的內(nèi)部，負(fù)載以列為單位進(jìn)行排序，要轉(zhuǎn)移的負(fù)載為從最大負(fù)載列開始，隔列轉(zhuǎn)移。這樣做的原因是由于燃燒的穩(wěn)定性使得相鄰網(wǎng)格列上的負(fù)載接近，通過排序并隔行轉(zhuǎn)移可以使得轉(zhuǎn)移的源進(jìn)程與目的進(jìn)程的負(fù)載相近，避免過多或過少轉(zhuǎn)移負(fù)載，引入新的負(fù)載不平衡。主控進(jìn)程標(biāo)記轉(zhuǎn)移的方向及數(shù)量，并廣播給所有計(jì)算進(jìn)程。然后根據(jù)上一步收到的標(biāo)記數(shù)組，負(fù)載小于平均負(fù)載的進(jìn)程（欠載進(jìn)程）準(zhǔn)備接收數(shù)組，并擴(kuò)大計(jì)算涉及的輸入及輸出數(shù)組，從過載進(jìn)程接收計(jì)算需要的輸入數(shù)組包，拆包將數(shù)據(jù)放入適當(dāng)?shù)奈恢?；過載進(jìn)程準(zhǔn)備發(fā)送數(shù)組，將計(jì)算需要的輸入數(shù)組打包，并發(fā)送給對應(yīng)的欠載進(jìn)程。

在計(jì)算完成之后，原來的欠載進(jìn)程將計(jì)算后改變的輸出數(shù)組打包，發(fā)送給原來的過載進(jìn)程。并恢復(fù)計(jì)算涉及的輸入及輸出數(shù)組為原來的規(guī)模。原來的過載進(jìn)程準(zhǔn)備接受數(shù)組包，并將其放入原先的位置。

5 效果的評價(jià)

本文從每一個(gè)計(jì)算步中，進(jìn)程的負(fù)載統(tǒng)計(jì)、有負(fù)載平衡策略與無負(fù)載平衡策略時(shí)程序的運(yùn)行時(shí)間、進(jìn)程之間的負(fù)載不平衡程度等方面對負(fù)載平衡的策略進(jìn)行分析。

（1）負(fù)載不平衡程度的改善

由于所有進(jìn)程中負(fù)載最大的進(jìn)程執(zhí)行時(shí)間最長，直接影響到程序執(zhí)行的總時(shí)間；而平均負(fù)載是理想情況下每個(gè)進(jìn)程的執(zhí)行時(shí)間。

有負(fù)載轉(zhuǎn)移和無負(fù)載轉(zhuǎn)移時(shí)，最大負(fù)載與平均負(fù)載的比值如圖9 所示，橫軸為程序執(zhí)行的步驟，灰色實(shí)線為無負(fù)載轉(zhuǎn)移時(shí)最大負(fù)載與平均負(fù)載的比值，虛線部分為使用牛頓迭代次數(shù)預(yù)測負(fù)載策略時(shí)最大負(fù)載與平均負(fù)載的比值，黑色實(shí)線部分為結(jié)合兩種預(yù)測方式后最大負(fù)載與平均負(fù)載的比值。從圖中可以看出，使用牛頓迭代次數(shù)預(yù)測在絕大部分時(shí)間可以較好地降低最大負(fù)載與平均負(fù)載的比值，但是在400 步左右由于化學(xué)反應(yīng)劇烈，這種策略就不夠有效。加入放熱量Q 預(yù)測負(fù)載后，在化學(xué)反應(yīng)劇烈的階段也能很好地預(yù)測負(fù)載變動(dòng)趨勢，最大負(fù)載與平均負(fù)載的比值降低到無負(fù)載轉(zhuǎn)移的一半左右。

圖9 最大負(fù)載與平均負(fù)載的比值

有負(fù)載轉(zhuǎn)移和無負(fù)載轉(zhuǎn)移時(shí)，各進(jìn)程的負(fù)載如圖10所示。

圖10 不同燃燒階段進(jìn)程上的負(fù)載對比

圖10的橫軸為各進(jìn)程編號，縱軸為該階段的累計(jì)負(fù)載，淺色直方圖為無負(fù)載轉(zhuǎn)移的統(tǒng)計(jì)負(fù)載，深色直方圖為有負(fù)載轉(zhuǎn)移的累計(jì)負(fù)載。圖10（a）顯示了劇烈化學(xué)反應(yīng)開始之前，在各進(jìn)程上無負(fù)載轉(zhuǎn)移和有負(fù)載轉(zhuǎn)移的負(fù)載，有負(fù)載轉(zhuǎn)移的部分負(fù)載更為平均。圖10（b）顯示了劇烈化學(xué)反應(yīng)開始之后，無負(fù)載轉(zhuǎn)移和有負(fù)載轉(zhuǎn)移的負(fù)載累計(jì)，有負(fù)載轉(zhuǎn)移的部分仍然更為平均。從圖10 還可以看出，有負(fù)載轉(zhuǎn)移的最大負(fù)載接近無負(fù)載轉(zhuǎn)移的最大負(fù)載的一半，進(jìn)程間的負(fù)載在有負(fù)載轉(zhuǎn)移時(shí)更加平衡。

（2）負(fù)載轉(zhuǎn)移時(shí)間分析

為了使負(fù)載轉(zhuǎn)移有效地執(zhí)行，即負(fù)載轉(zhuǎn)移的開銷不能大于負(fù)載轉(zhuǎn)移帶來的性能提升，因此需要度量負(fù)載轉(zhuǎn)移帶來的開銷。

由于每一個(gè)迭代步中，影響程序執(zhí)行時(shí)間的部分為負(fù)載最大的進(jìn)程的執(zhí)行之間，因此使用每個(gè)迭代步負(fù)載最大的進(jìn)程中，負(fù)載轉(zhuǎn)移開銷占總的求解常微分方程計(jì)算時(shí)間的比值作為負(fù)載轉(zhuǎn)移的開銷。如圖11 所示。

圖11 負(fù)載轉(zhuǎn)移所占時(shí)間比例

由圖11 中可以看出，負(fù)載轉(zhuǎn)移所占時(shí)間比例最大不超過總執(zhí)行時(shí)間的35%，而根據(jù)負(fù)載轉(zhuǎn)移策略，是將過載進(jìn)程與欠載進(jìn)程負(fù)載差值的一半轉(zhuǎn)移出去，即計(jì)算時(shí)間的50%，即節(jié)省了50%的時(shí)間，而最大的開銷僅為30%，特別是劇烈化學(xué)反應(yīng)負(fù)載比較大的情況下，負(fù)載轉(zhuǎn)移策略的開銷不到計(jì)算時(shí)間的5%，由此可以得出負(fù)載轉(zhuǎn)移一定是有效的。

（3）程序總體性能的提升

為了分析負(fù)載轉(zhuǎn)移策略對主程序帶來的性能提升，需要分析求解常微分方程計(jì)算部分占整個(gè)程序運(yùn)行時(shí)間的百分比。經(jīng)過實(shí)驗(yàn)得出，在程序的核心計(jì)算過程中，解常微分方程計(jì)算時(shí)間占總執(zhí)行時(shí)間的5%～35%。

定義性能提升的度量為在并行情況下，不同核數(shù)執(zhí)行程序核心計(jì)算過程時(shí)，和有負(fù)載平衡的執(zhí)行時(shí)間與無負(fù)載平衡的執(zhí)行時(shí)間減少時(shí)間的比值。測試網(wǎng)格初始時(shí)為24×4，化學(xué)反應(yīng)穩(wěn)定時(shí)的規(guī)模在30 000 左右。

從表1中可以看出，程序的核心計(jì)算性能總體隨著核數(shù)的增加而提高，12 核運(yùn)行時(shí)程序性能提升達(dá)到了10.20%。

表1 程序的性能提升

6 總結(jié)和未來工作

本文對非結(jié)構(gòu)網(wǎng)格自適應(yīng)的燃燒并行數(shù)值模擬計(jì)算的特點(diǎn)進(jìn)行了詳細(xì)的研究，提出了數(shù)據(jù)按列存儲的改進(jìn)方式和根據(jù)不同化學(xué)反應(yīng)情況預(yù)測負(fù)載的機(jī)制，還提出了將負(fù)載一對一地轉(zhuǎn)移到欠載進(jìn)程上動(dòng)態(tài)負(fù)載平衡策略，將最大負(fù)載與平均負(fù)載的比值降低了一半左右，程序總體性能提升在9.00%～10.20%。由于一對一負(fù)載轉(zhuǎn)移的限制，還有部分進(jìn)程處于空閑狀態(tài)，下一步的工作是將負(fù)載轉(zhuǎn)移改為一對多，將過載進(jìn)程上的負(fù)載盡可能多地分出去，使得所有進(jìn)程上的負(fù)載更加平衡。

[1] Chen Z，Burke M P，Ju Y.Effects of compression and stretch on the determination of laminar flame speeds using propagating spherical flames[J].Combustion Theory and Modelling，2009，13（2）：343-364.

[2] Sun，Takayama.Conservative smoothing on an adaptive quadrilateral grid[J].Journal of Computational Physics，1999，150（1）：143-180.

[3] Hummel S F，Schonberg E，F(xiàn)lynn L E.Factoring：a method for scheduling parallel loops[J].Comm of the ACM，1992，35（8）：90-101.

[4] Kumar V，Grama A Y，RaoVempaty N.Scalable load balancing techniques for parallel computers[J].J of Parallel and Distr，1994，22（1）：60-79.

[5] Lain M S，Rothberg E E，Wolf M E.The cache performance and optimizations of block algorithms[C]//4th Int Conf on Architectural Support for Progr Lang and Operating Systems，1991：63-74.

[6] Orlando S，Perego R.A template for non-uniform parallel loops based on dynamic scheduling and prefetching techniques[C]//Proc of the 1996 ACM Int Conf on Supercomputing，1996：117-124.

[7] Orlando S，Perego H.SUPPLE：an efficient run-time support for non-uniform parallel loops[J].J of System Architecture，1999，45（15）：1323-1343.

[8] Oliker L，Biswas R.PLUM：Parallel Load balancing for adaptive Unstructured Meshes[J].Parallel and Distributed Computing，1998，52（2）：150-177.

[9] Bui T.Graph bisection algorithms with good average behavior[J].Combinatorica，1984，7（2）：181-191.

[10] Kernighan B，Lin S.An efficient heuristic procedure for partitioning graphs[J].Bell System Technical Journal，1970，49（2）：291-307.

[11] Chan T，Mathew T.Domain decomposition algorithms[J].Acta Numerica，1994，3（1）：61-143.

[12] 王學(xué)德，伍貽兆，夏健.動(dòng)態(tài)負(fù)載平衡的二維非結(jié)構(gòu)網(wǎng)格DSMC并行算法研究[J].空氣動(dòng)力學(xué)學(xué)報(bào)，2007，25（3）：339-361.

[13] 劉鑫，陸林生，陳德訓(xùn).非結(jié)構(gòu)網(wǎng)格并行計(jì)算預(yù)處理方法研究[J].計(jì)算機(jī)科學(xué)，2012，39（3）：308-311.

[14] 陳建軍.非結(jié)構(gòu)網(wǎng)格的并行劃分和自適應(yīng)非結(jié)構(gòu)化網(wǎng)格生成及其并行化的若干問題研究[D].杭州：浙江大學(xué)，2006.

[15] 李桂波，楊國偉.基于多塊結(jié)構(gòu)網(wǎng)格的并行計(jì)算及負(fù)載平衡研究[J].宇航學(xué)報(bào)，2011，32（6）：1224-1230.