高效稀疏矩陣存儲及油藏數值模擬的方法和裝置的制造方法
【專利摘要】本發明公開了一種高效稀疏矩陣存儲及油藏數值模擬的方法和裝置，屬于油藏數值模擬技術領域，本發明采用新型稀疏矩陣存儲格式，其典型特征是對角元與非對角元分開存儲，元素首先按塊排列，塊內按列排列，下標數組利用已有網格信息；同時提供對應的“矩陣?向量相乘”算法。本發明中，新型稀疏矩陣存儲格式有效利用了油藏矩陣的塊狀特性、非零元位置對稱特性以及非對角元位置等同于連通表的特性，使得矩陣存儲空間有效減小；由于存儲格式的特點，“矩陣?向量相乘”算法可以有效利用計算機緩存，加速計算的進行，最終對加速整個數值模擬過程。
【專利說明】
高效稀疏矩陣存儲及油藏數值模擬的方法和裝置
技術領域
[0001]本發明涉及油藏數值模擬技術領域，特別是指一種高效稀疏矩陣存儲及油藏數值模擬的方法和裝置。
【背景技術】
[0002]油藏數值模擬技術是評價油藏潛力，制定開發方案的基礎。并且在現代油藏開發中已經有了廣泛的使用。油藏數值模擬的過程就是將物理問題用數學方程刻畫，經過一系列的離散求解，最終得到想要的結果。
[0003]在求解過程的最后，需要求解一個大型線性方程組，此方程組可以用矩陣形式表示，形如Ax = b，這里A就代表矩陣，一般是一個大型稀疏塊狀矩陣(如圖1所示，黑點代表非零元素)。對于這種大型稀疏矩陣，有兩個重要的問題:一是如何存儲;一是相應的“矩陣-向量相乘”操作是否高效。
[0004]大型稀疏矩陣在很多物理模擬問題中都會出現，其特點是矩陣總元素數極其龐大(可能大于18個)，而其中的非零元素數只占很小的一個比例。通用的稀疏矩陣存儲方法都是只存儲其中的非零元素，再以不同的方式存儲這些非零元素在原矩陣中的位置信息，典型的存儲格式有壓縮行存儲(CSR)、壓縮列存儲(CSC)、坐標存儲(COO)等。
[0005]“矩陣-向量相乘”是求解大型線性方程組的迭代算法中經常使用到的一個操作步驟。這一步操作能否高效完成直接影響到整個求解算法的效率，而這步操作如何進行以及是否快捷又與具體選擇的矩陣存儲格式息息相關。對于主流的通用稀疏矩陣存儲格式，網上很多矩陣運算庫中都提供了非常優化的“矩陣-向量相乘”函數(如SparseLib等)。
[0006]在油藏數值模擬中，可以使用這些通用的矩陣存儲格式以及相應的“矩陣-向量相乘”算法，但它們的問題在于它們并非針對油藏數值模擬而設計，因此沒有考察油藏矩陣的特殊性質，不能達到最高的效率。
【發明內容】
[0007]本發明要解決的技術問題是提供一種高效稀疏矩陣存儲及油藏數值模擬的方法和裝置，其能有效節省矩陣存儲空間，加快矩陣生成速度，提供高效的“矩陣-向量相乘”運算，從而加速整個油藏數值模擬的進行。
[0008]為解決上述技術問題，本發明提供技術方案如下:
[0009]—方面，提供一種高效稀疏矩陣存儲的方法，包括:
[0010]數值數組存儲:I)將矩陣非零元素根據正方塊狀對角元和上下非對角元分別存儲在三個數組中，分別用mDiag，m0ffDiag_B表示；2)每個數組內部，按塊順序排列，對角元數組mDiag按照塊從上到下排列，上下非對角元數組應位置元素塊也按對稱方式存儲;對mOffDiag_A而言，按照從上到下，從左到右的順序排列，也就是按照和對角元塊平行的方式，由近至遠排列;1110打01&8_8每個相應位置的塊都和mOffDiag_A*的塊對稱;3)每個塊內按列優先的順序排列，即在每一塊內，按照從上到下，從左到右的順序排列；
[0011 ]下標數組存儲:無需額外存儲下標數組。
[0012]另一方面，提供一種基于上述的高效稀疏矩陣存儲方法的油藏數值模擬的方法，包括:
[0013]步驟1:選取一種稠密矩陣向量相乘算法；
[0014]步驟2:按照存儲順序，依次做對角元數組的乘法；
[0015]步驟3:按照存儲順序，依次做兩個非對角元數組的乘法，并且累加對應的對角元數組的乘法結果，即得最終的結果。
[0016]進一步的，所述稠密矩陣向量相乘算法為InteI的MKL算法。
[0017]再一方面，提供一種高效稀疏矩陣存儲的裝置，包括存儲模塊，用于:
[0018]數值數組存儲:I)將矩陣非零元素根據正方塊狀對角元和上下非對角元分別存儲在三個數組中，分別用mDiag，m0ffDiag_B表示；2)每個數組內部，按塊順序排列，對角元數組mDiag按照塊從上到下排列，上下非對角元數組應位置元素塊也按對稱方式存儲;對mOffDiag_A而言，按照從上到下，從左到右的順序排列，也就是按照和對角元塊平行的方式，由近至遠排列;1110打01&8_8每個相應位置的塊都和mOffDiag_A*的塊對稱;3)每個塊內按列優先的順序排列，即在每一塊內，按照從上到下，從左到右的順序排列；
[0019]下標數組存儲:無需額外存儲下標數組。
[0020]又一方面，提供一種基于上述的高效稀疏矩陣存儲裝置的油藏數值模擬的裝置，包括:
[0021]選取模塊，用于選取一種稠密矩陣向量相乘算法；
[0022]第一乘法模塊，用于按照存儲順序，依次做對角元數組的乘法；
[0023]第二乘法模塊，用于按照存儲順序，依次做兩個非對角元數組的乘法，并且累加對應的對角元數組的乘法結果，即得最終的結果。
[0024]進一步的，所述稠密矩陣向量相乘算法為InteI的MKL算法。
[0025]本發明具有以下有益效果:
[0026]上述方案中，定義了新型稀疏矩陣存儲格式(后文稱作CIB格式)，其典型特征是對角元與非對角元分開存儲，元素首先按塊排列，塊內按列排列，下標數組利用已有網格信息;提供對應的“矩陣-向量相乘”算法(后文稱作SpMV算法)。
[0027]本發明中，CIB存儲格式有效利用了油藏矩陣的塊狀特性、非零元位置對稱特性以及非對角元位置等同于連通表的特性，使得矩陣存儲空間有效減小；由于存儲格式的特點，基于CIB的SpMV算法可以有效利用計算機緩存，加速計算的進行，最終對加速整個數值模擬過程。
【附圖說明】
[0028]圖1為典型的油藏數值模擬問題產生的大型稀疏矩陣結構；
[0029]圖2為簡單的包含6個網格的油藏系統，以及其對應連通表；
[0030]圖3為本發明的高效稀疏矩陣存儲格式(CIB格式)具體的存儲方法示意圖；
[0031]圖4為本發明的“矩陣-向量相乘”算法(SpMV算法)示意圖；
[0032]圖5為不同矩陣塊(block)大小及網格規模下CIB格式對CSR格式SpMV運算的加速比(CSR耗時除以CIB耗時)。
【具體實施方式】
[0033]為使本發明要解決的技術問題、技術方案和優點更加清楚，下面將結合附圖及具體實施例進行詳細描述。
[0034]本申請的發明人經過潛心研究發現:在油藏數值模擬系統中，獲得的大型稀疏矩陣的一個顯著特征是其呈現出塊狀分布:每一個小塊都與一個數值模擬網格相對應，小塊的行數代表對應網格擁有的方程個數，小塊的列數代表對應網格擁有的變量個數。對常見的全隱式解法而言，方程數等于變量數，也就是說每一個小塊都是行列相等的正方塊。
[0035]整個矩陣本質上是一個偏導矩陣。塊狀對角元部分代表網格i的方程對網格i自身變量的偏導，所以每一個塊狀對角元總是存在的;塊狀非對角元部分代表網格i的方程對網格j變量的偏導，所以如果網格i和網格j是相鄰的，那么在塊狀矩陣的(i，j)位置和其對稱位(j，i)位置的矩陣塊就存在，反之，如果i，j不相鄰，這兩個位置就是O。哪些網格是相連的，這一信息在數值模擬程序中一般由網格處理模塊分析地質網格模型后導出，存儲在“連通表”(又稱“鏈接表”)中。
[0036]任何一種稀疏矩陣存儲方法都要提供數值數組和下標數組。數值數組的本質是壓縮一一只存儲非零元而剔除零元;數組壓縮后喪失了原矩陣的位置信息，這時就需要下標數組來告知用戶壓縮后數組中的每一個元素在原矩陣中是第幾行第幾列。
[0037]下面將詳細說明本發明的高效稀疏矩陣存儲方法的實現方式。首先介紹數值數組如何存儲，然后是下標數組，最后是此方法的好處。
[0038]數值數組存儲:
[0039]I)將矩陣非零元素根據正方塊狀對角元和上下非對角元分別存儲在三個數組中，分別用mDiag表示，其中mDiag代表塊狀對角元存入的數組，
分別代表對角元上部的非對角元和對角元下部的非對角元存入的數組。
[0040]2)每個數組內部，首先按塊順序排列。對角元數組mDiag就按照塊從上到下排列。上下非對角元是對稱的，所以對應位置元素塊也按對稱方式存儲。對mOffDiag_A而言，按照從上到下，從左到右的順序排列，也就是按照和對角元塊平行的方式，由近至遠排列c^mOf個相應位置的塊都和mOffDiag_A*的塊對稱。
[0041 ] 3)每個塊內按列優先的順序排列，即在每一塊內，按照從上到下，從左到右的順序排列。至此矩陣非零元素的存儲排列規則已完全說明。
[0042]下標數組存儲:
[0043]在我們的這種方法中，由于可以重復利用已有信息，下標數組可以節省，無需額外存儲。因為是按塊狀存儲，只需要確定每一個塊的位置，再加上已知的塊大小，就可以推得每一個塊內元素在原矩陣中的位置。對于對角元，由于對角元數組按從上到下的順序排列，它在數組中的位置就是它在原矩陣中的行和列；對于非對角元，mOf
數組相同位置的元素在原矩陣中是對稱的，只需給出一個的下標，交換行號和列號就可得到另一個的下標，非對角元數組的排列順序是經過特別選取的，這個順序和連通表中記錄的網格相連的信息順序一致，可以直接使用連通表來得知非對角元數組中元素在原矩陣中的位置，無須額外存儲。
[0044]新存儲方法的好處:I)從存儲效率上而言，利用了油藏問題產生的矩陣的塊狀結構和對稱性，使得可以利用已有信息使下標數組的存儲大大減少。2)從矩陣生成效率而言，傳統存儲模式比較抽象，壓縮過后的矩陣物理含義不明確，所以不容易直接生成，需要一個額外的“壓縮”步驟，即計算出的各網格方程對各變量偏導后先存儲為易于理解的形式，再對此“矩陣”進行壓縮得到最后的壓縮存儲矩陣。而我們提出的新的存儲模式每一個元素都有明確的物理意義，可以與計算偏導的過程直接對應，直接生成壓縮矩陣形式。3)計算效率而言，新存儲模式的特征允許我們為其設計更高效的“矩陣一向量相乘”算法。
[0045]同時，本發明提供一種基于上述的高效稀疏矩陣存儲方法的油藏數值模擬的方法，即對應此新型存儲格式的“矩陣-向量相乘”算法，包括:
[0046]步驟1:選取一種稠密矩陣向量相乘算法，這個算法用于后面矩陣小塊與向量的乘法計算。此類算法已經得到充分的研究，可以從很多矩陣運算庫中獲取優化過的函數，比如Intel的MKL;
[0047]步驟2:按照存儲順序，依次做對角元數組的乘法。
[0048]本步驟中，假設mDiag代表對角元數組，mDiag_i代表對角元數組中第i個矩陣塊，V代表要相乘的向量，V_i代表對應大小的第i個向量塊，W代表結果向量，W_i代表對應大小的第i個向量塊。那么對角元數組的乘法可以表示為:按照順序，對每一個i，ff_i = mDiag_i XV」。
[0049]步驟3:按照存儲順序，依次做兩個非對角元數組的乘法，并且累加對應的對角元數組的乘法結果，即得最終的結果。
[0050]本步驟中，上下非對角元處理方法相似，只以上對角元為例。假設H1ffDiagAjR表上非對角元中第i個矩陣塊，(ai，bi)代表連通表中存儲的第i對鏈接，V_i和W_i意義同步驟2。那么上非對角元數組的乘法可以表示為:對于每一個i，W_ai + = mOffDiagA_i X V_bi(“+=”符號表示在之前的結果上累加，因為W中已經有對角元數組乘法的結果，需要累加非對角元的乘法結果才是最終的結果)。
[0051]下面將結合一個具體的實例予以說明。現有一個6個網格(呈2行3列分布)的油藏系統，其示意圖以及對應的連通表如圖2所示。對于兩相流模擬，從這個問題導出的矩陣如圖3中第一幅圖所示。整個矩陣由多個2X2的小矩陣塊組成。
[0052]圖中11^&8代表塊狀對角元存入的數組，111(^0丨&8_4和111(^0丨&8_8分別代表對角元上部的非對角元和對角元下部的非對角元存入的數組。圖中對角元數組mD i a g就首先按照塊從上到下排列，即先是(0，0)塊，然后是(I，1)塊，以此類推。上下非對角元是對稱的，所以111(^€0丨&8_4和111(^€0丨&8_8對應位置元素塊也按對稱方式存儲。對111(^€0丨&8_4而言，按照從上到下，從左到右的順序排列，也就是按照和對角元塊平行的方式，由近至遠排列。先是(0，1)塊，然后是(1，2)塊，然后是(3，4)塊，以此類推。對稱地，1110打01&8_8先存(1，0)塊，然后是(2，I)塊，然后是(4，3)塊，等等。如圖中帶箭頭的曲線所示，在每一塊內，按照從上到下，從左到右的列優先順序排列。
[0053]對于mDiag數組，不需要下標數組輔助，即可得知每一個元素在原矩陣中的位置。以其中第7個元素為例。首先算得它所在矩陣塊在原矩陣中的位置，7除以4等于I余3，所以它位于第2個塊，在第2個塊之前一共有1*2 = 2行，1*2 = 2列。它在塊內的第I行第2列，所以綜合起來，它在原矩陣中是第2+1 = 3行，第2+2 = 4列。
[0054]對于mOf fDiag_A數組，連通表中的信息等同于其位置信息，所以下標數組也可以省略。以其中第21個元素為例。首先算得它所在矩陣塊在原矩陣中的位置，21除以4等于5余I，所以它位于第6個塊。在連通表中查找第6組連接信息，為I，4，所以它所在的矩陣塊在原塊狀矩陣的第2行第5列(注意編號從O開始)。在這個塊之前一共有1*2 = 2行，4*2 = 8列。它在塊內的第I行第I列，所以綜合起來，它在原矩陣中是第2+1 = 3行，第8+1=9列。m0ffDiag_B中塊與mOf f Diag_A對稱，同理可得。
[0055]下面繼續結合例子說明“矩陣-向量相乘”算法的實施過程。假設A是圖4中矩陣，X是要乘的向量，b是乘得的結果向量。圖中各個矩陣塊已標出。計算過程為:
[0056]首先，對角元數組與向量b相乘。按順序，對每一個對角元塊Di(i= 0,1，2，…5)，
[0057]bi=Di Xxi
[0058]即b0= D0 X x0，bl = Dl X xl，b2 = D2 X x2等等。
[0059]然后，上非對角元數組與向量b相乘。按連通表順序，對每一個上非對角元塊OAi(i=1，2，...7)，找到連通表中相應連接信息(r，c)，計算
[0060]br+ = 0Ai X xc
[0061 ]即b0+ = 0Al X xl，bl+ = 0A2 X x2，b3+ = 0A3 X x4,等等。
[0062]最后，下非對角元數組與向量b相乘，由于對稱性，按上非對角元做法，交換行列，即r和c即可。
[0063]基于新存儲結構的“矩陣-向量相乘”算法利用了油藏問題矩陣的特征，其算法也可以更有效地利用計算機緩存，因此相比傳統存儲格式下的算法速度上有很大提高。我們的測試結果見圖5，對大部分算例來說加速比在2倍以上。
[0064]另一方面，與上述的方法相對應，本發明還提供一種高效稀疏矩陣存儲的裝置，包括存儲模塊，用于:
[0065]數值數組存儲:I)將矩陣非零元素根據正方塊狀對角元和上下非對角元分別存儲在三個數組中，分別用mDiag，m0ffDiag_B表示；2)每個數組內部，按塊順序排列，對角元數組mDiag按照塊從上到下排列，上下非對角元數組應位置元素塊也按對稱方式存儲;對mOffDiag_A而言，按照從上到下，從左到右的順序排列，也就是按照和對角元塊平行的方式，由近至遠排列;1110打01&8_8每個相應位置的塊都和mOffDiag_A*的塊對稱;3)每個塊內按列優先的順序排列，即在每一塊內，按照從上到下，從左到右的順序排列；
[0066]下標數組存儲:無需額外存儲下標數組。
[0067]同時，還提供一種基于上述的高效稀疏矩陣存儲裝置的油藏數值模擬的裝置，包括:
[0068]選取模塊，用于選取一種稠密矩陣向量相乘算法；
[0069]第一乘法模塊，用于按照存儲順序，依次做對角元數組的乘法；
[0070]第二乘法模塊，用于按照存儲順序，依次做兩個非對角元數組的乘法，并且累加對應的對角元數組的乘法結果，即得最終的結果。
[0071]進一步的，所述稠密矩陣向量相乘算法為InteI的MKL算法。
[0072]以上所述是本發明的優選實施方式，應當指出，對于本技術領域的普通技術人員來說，在不脫離本發明所述原理的前提下，還可以作出若干改進和潤飾，這些改進和潤飾也應視為本發明的保護范圍。
【主權項】
1.一種高效稀疏矩陣存儲的方法，其特征在于，包括: 數值數組存儲:I)將矩陣非零元素根據正方塊狀對角元和上下非對角元分別存儲在三個數組中，分別用mDiag，mOffDiag_B表示；2)每個數組內部，按塊順序排列，對角元數組mDiag按照塊從上到下排列，上下非對角元數組對應位置元素塊也按對稱方式存儲;對mOffDiag_A而言，按照從上到下，從左到右的順序排列，也就是按照和對角元塊平行的方式，由近至遠排列；個相應位置的塊都和mOffDiag_A*的塊對稱;3)每個塊內按列優先的順序排列，即在每一塊內，按照從上到下，從左到右的順序排列； 下標數組存儲:無需額外存儲下標數組。2.基于權利要求1所述的高效稀疏矩陣存儲方法的油藏數值模擬的方法，其特征在于，包括: 步驟1:選取一種稠密矩陣向量相乘算法； 步驟2:按照存儲順序，依次做對角元數組的乘法； 步驟3:按照存儲順序，依次做兩個非對角元數組的乘法，并且累加對應的對角元數組的乘法結果，即得最終的結果。3.根據權利要求2所述的油藏數值模擬的方法，其特征在于，所述稠密矩陣向量相乘算法為Inte I的MKL算法。4.一種高效稀疏矩陣存儲的裝置，其特征在于，包括存儲模塊，用于: 數值數組存儲:I)將矩陣非零元素根據正方塊狀對角元和上下非對角元分別存儲在三個數組中，分別用mDiag，mOffDiag_B表示；2)每個數組內部，按塊順序排列，對角元數組mDiag按照塊從上到下排列，上下非對角元數組對應位置元素塊也按對稱方式存儲;對mOffDiag_A而言，按照從上到下，從左到右的順序排列，也就是按照和對角元塊平行的方式，由近至遠排列；個相應位置的塊都和mOffDiag_A*的塊對稱;3)每個塊內按列優先的順序排列，即在每一塊內，按照從上到下，從左到右的順序排列； 下標數組存儲:無需額外存儲下標數組。5.基于權利要求4所述的高效稀疏矩陣存儲裝置的油藏數值模擬的裝置，其特征在于，包括: 選取模塊，用于選取一種稠密矩陣向量相乘算法； 第一乘法模塊，用于按照存儲順序，依次做對角元數組的乘法； 第二乘法模塊，用于按照存儲順序，依次做兩個非對角元數組的乘法，并且累加對應的對角元數組的乘法結果，即得最終的結果。6.根據權利要求5所述的油藏數值模擬的裝置，其特征在于，所述稠密矩陣向量相乘算法為Inte I的MKL算法。
【文檔編號】G06F17/50GK105844009SQ201610165195
【公開日】2016年8月10日
【申請日】2016年3月22日
【發明人】龔斌, 石欣, 張忠國
【申請人】北京大學