復雜介質彈性波傳播模擬的多網(wǎng)格切比雪夫并行譜元法
【專利摘要】本發(fā)明公開了一種多網(wǎng)格切比雪夫并行譜元方法。該方法包括：a.將計算區(qū)域劃分為大而規(guī)則的單元；b.在單元內(nèi)定義主網(wǎng)格和輔助網(wǎng)格，主網(wǎng)格為切比雪夫配置點；c.利用主網(wǎng)格將單元內(nèi)波場作切比雪夫截斷展開，利用輔助網(wǎng)格將單元內(nèi)介質參數(shù)和外力作截斷展開；d.將波場、介質參數(shù)和外力的截斷展開式代入波動方程，得到單元質量矩陣、單元剛度矩陣和單元外力向量；e.基于各單元的單元質量矩陣、單元剛度矩陣和單元外力向量，在主網(wǎng)格上求解彈性波方程。
【專利說明】復雜介質彈性波傳播模擬的多網(wǎng)格切比雪夫并行譜元法
【技術領域】
[0001]本發(fā)明涉及地球物理，工程地震學，計算聲學等領域的高性能計算。
【背景技術】
[0002]復雜介質中彈性波的傳播模擬在地球物理，工程地震學和計算聲學等領域具有重要地位。如何找到一種更精確、更快速、更適宜并行化的數(shù)值模擬方法也一直是國、內(nèi)外相關領域研究者的工作目標和重點。
[0003]當前在彈性波方程模擬中普遍使用的數(shù)值方法，主要有有限差分法、有限元法、偽譜法和譜元法等。它們有各自的優(yōu)缺點并在不同的領域得到應用。
[0004]有限差分法將方程中的導數(shù)用差分近似，從而將偏微分方程轉換為代數(shù)方程來求解，其誤差由離散化配置點數(shù)目和Taylor級數(shù)的截斷誤差決定。有限差分法由于編程簡捷和計算速度較快被廣泛應用于計算地球物理和工程地震中。但低階有限差分法需要大量的網(wǎng)格點來保證精度，而高階有限差分法又很難高效地處理自由界面和復雜結構問題。
[0005]有限元方法基于波動方程的弱形式，將計算區(qū)域劃分為有限個互不重疊的單元，在每個單元內(nèi)，波場由低階多項式(如分段線性函數(shù))來近似。有限元法適合處理邊界和內(nèi)部復雜結構，在結構分析和瞬態(tài)模擬方面獲得了廣泛的應用。但有限元法在大型的彈性波方程模擬領域應用相對較少，其原因一是有限元法需要耗用大量計算資源，二是低階有限元存在頻散現(xiàn)象，而通常的高階有限元存在偽波。
[0006]偽譜法利用無限可微的全局函數(shù)(如富里葉Fourier級數(shù))，將問題從物理空間轉化為波數(shù)域進行處理，具有良好的處斂性。偽譜法只需要很少的空間網(wǎng)格點數(shù)即可達到很高的精度。為了處理邊界條件，在空間展開時經(jīng)常使用切比雪夫Chebyshev或拉格朗日Legendre正交多項式來取代Fourier級數(shù)。但由于正交多項式的配置點分布很不均勻，給時間步長的選擇帶來很大限制。同時，如同有限差分法一樣，偽譜法不能自如地處理彎曲的自由界面和復雜結構的問題。為解決這類問題有人提出了彎曲坐標和區(qū)域分解等方法，但代價是計算量的增加。
[0007]譜元法由Patera在1984年提出，早期主要應用于流體力學。Seriani等在1991年首次將Chebychev譜元法引入到聲波方程的數(shù)值模擬中，各國研究者在其后做了大量的相關研究。它把有限元法和偽譜法相結合，兼具了有限元處理邊界和結構的靈活性和偽譜法的快速收斂特性。在達到相同的精度前提下，可以采用比傳統(tǒng)有限元更稀疏的單元劃分，減少了計算時間和內(nèi)存需求。
[0008]譜元法的基本方法是在每一個單元上用高階譜展開。選取以截斷的正交多項式表示的基函數(shù)，在各個單元上利用配置點插值，以提高解的收斂速度。其主要步驟是:(1)首先把計算的區(qū)域分成許多子域(單元)，每個子域由若干節(jié)點(配置點)組成；(2)在每個子域中把近似解表示成截斷的正交多項式展開；(3)用Galerkin方法求解正交問題的變分格式，得到全局的近似解。
[0009]譜元法可以用較稀疏的網(wǎng)格和單元獲得較高精度，但傳統(tǒng)的譜元法中，每個單元內(nèi)只能有單一均勻介質，某些情況下嚴重降低計算效率。比如當介質結構復雜，變化尺度很小甚至小于波長時，必須按介質變化的小尺度劃分單元和求解，這樣將會造成極大的計算規(guī)模。
[0010]此外，傳統(tǒng)的譜元法求解波動方程時，需要形成全局剛度矩陣和質量矩陣，考慮衰減時還需引進阻尼矩陣，需要耗費大量的存貯空間，在算法上限制了其并行化的效率。注意到在波動方程的時間迭代過程中，需要做多次矩陣和向量乘積，如Spk，這里S和Pk分別為全局剛度(或質量)矩陣和全局向量。考慮到全局矩陣的稀疏特性，全局的矩陣和相量相乘耗費了大量的計算時間。
[0011]以上兩點都極大地阻礙了譜元法在地球物理，工程地震學，計算聲學等領域的大規(guī)模實際應用。

【發(fā)明內(nèi)容】

[0012]本發(fā)明的目的在于提供一種能夠克服上述缺點的算法。
[0013]為此，本發(fā)明提供一種多網(wǎng)格Chebyshev并行譜元方法。該方法包括:
[0014]a.將計算區(qū)域劃分為大而規(guī)則的單元；
[0015]b.在單元內(nèi)定義主網(wǎng)格和輔助網(wǎng)格，主網(wǎng)格為Chebyshev配置點；
[0016]c.利用主網(wǎng)格將單元內(nèi)波場作Chebyshev截斷展開，利用輔助網(wǎng)格將單元內(nèi)介質參數(shù)和外力作截斷展開；
[0017]d.將波場、介質參數(shù)和外力的截斷展開式代入波動方程，得到單元質量矩陣、單元剛度矩陣和單元外力向量；
[0018]e.基于各單元的單元質量矩陣、單元剛度矩陣和單元外力向量，在主網(wǎng)格上求解彈性波方程。
[0019]優(yōu)選地，利用預條件共軛梯度法求解彈性波方程。
[0020]進一步優(yōu)選地，彈性波方程包括全局矩陣與全局向量的乘積；所述利用預條件共軛梯度法求解彈性波方程的步驟包括:計算單元矩陣，從全局向量中按全局節(jié)點與單元節(jié)點的對應關系取出單元向量；在單元上獨立計算單元矩陣向量乘積，再將結果的單元向量按單元節(jié)點與全局節(jié)點的對應關系疊加，形成全局結果向量。
[0021 ] 進一步優(yōu)選地，將所有單元均分給P個CPU計算節(jié)點，所述利用預條件共軛梯度法求解彈性波方程的步驟包括:每個計算節(jié)點的進程獨立讀入該區(qū)域的介質參數(shù)，計算單元矩陣；通過疊加和傳遞子區(qū)域之間相鄰單元共有節(jié)點上的向量值，獲得全局向量在計算節(jié)點上分配的部分；在計算節(jié)點內(nèi)部計算單元矩陣和向量積，得到單元結果向量；疊加和傳遞相鄰單元共有節(jié)點上的結果向量值，得到全局結果向量在節(jié)點上分配的部分。
[0022]優(yōu)選的，基于各單元的單元質量矩陣、單元剛度矩陣和單元外力向量，在主網(wǎng)格上求解彈性波方程的步驟包括:在各進程中根據(jù)單元質量矩陣、單元剛度矩陣獨立計算矩
陣
【權利要求】
1.一種多網(wǎng)格切比雪夫Chebyshev并行譜元方法，其特征包括:a.將計算區(qū)域劃分為大而規(guī)則的單元 ；b.在單元內(nèi)定義主網(wǎng)格和輔助網(wǎng)格，主網(wǎng)格為Chebyshev配置點；c.利用主網(wǎng)格將單元內(nèi)波場作Chebyshev截斷展開，利用輔助網(wǎng)格將單元內(nèi)介質參數(shù)和外力作截斷展開；d.將波場、介質參數(shù)和外力的截斷展開式代入波動方程，得到單元質量矩陣、單元剛度矩陣和單元外力向量；e.基于各單元的單元質量矩陣、單元剛度矩陣和單元外力向量，在主網(wǎng)格上求解彈性波方程。
2.如權利要求1所述的多網(wǎng)格Chebyshev并行譜元方法，其中利用預條件共軛梯度法求解彈性波方程。
3.如權利要求2所述的多網(wǎng)格Chebyshev并行譜元方法，其中彈性波方程包括全局矩陣與全局向量的乘積；所述利用預條件共軛梯度法求解彈性波方程的步驟包括:計算單元矩陣，從全局向量中按全局節(jié)點與單元節(jié)點的對應關系取出單元向量；在單元上獨立計算單元矩陣向量乘積，再將結果的單元向量按單元節(jié)點與全局節(jié)點的對應關系疊加，形成全局結果向量。
4.如權利要求2所述的多網(wǎng)格Chebyshev并行譜元方法，其中，將所有單元均分給P個CPU計算節(jié)點，所述利用預條件共軛梯度法求解彈性波方程的步驟包括:每個計算節(jié)點的進程獨立讀入該區(qū)域的介質參數(shù)，計算單元矩陣；通過疊加和傳遞子區(qū)域之間相鄰單元共有節(jié)點上的向量值，獲得全局向量在計算節(jié)點上分配的部分；在計算節(jié)點內(nèi)部計算單元矩陣和向量積，得到單元結果向量；疊加和傳遞相鄰單元共有節(jié)點上的結果向量值，得到全局結果向量在節(jié)點上分配的部分。
5.如權利要求1所述的多網(wǎng)格Chebyshev并行譜元方法，其中，基于各單元的單元質量矩陣、單元剛度矩陣和單元外力向量，在主網(wǎng)格上求解彈性波方程的步驟包括:在各進程中根據(jù)單元質量矩陣、單元剛度矩陣獨立計算矩陣
6.如權利要求5所述的多網(wǎng)格Chebyshev并行譜元方法，其中，用預條件共軛梯度法結合逐元技術的步驟包括:各進程獨立初始化解，計算殘余向量；各進程傳遞相鄰單元共有節(jié)點對應的殘余向量值并疊加，得到更新的殘余向量；各進程獨立計算矩陣向量乘積，用并行算法計算向量內(nèi)積；各進程傳遞相鄰單元共有節(jié)點上的向量值并疊加，得到更新的矩陣向量乘積；更新解和殘余向量。
7.如權利要求1所述的多網(wǎng)格Chebyshev并行譜元方法，其中，每個單元內(nèi)部可含有多種介質。`
【文檔編號】G06F17/50GK103530451SQ201310452393
【公開日】2014年1月22日 申請日期:2013年9月27日 優(yōu)先權日:2013年9月27日
【發(fā)明者】蘇暢, G·塞里安尼, 林偉軍 申請人:中國科學院聲學研究所