基于圖卷積的一維平流方程空間離散化數值求解加速方法?

時津津, 宋 寧, 田 浩, 聶 婕??, 魏志強

(1.中國海洋大學信息科學與工程學部, 山東 青島 266100; 2.中國海洋大學數學科學學院, 山東 青島 266100)

隨著計算機技術的發展,偏微分方程在物理學、氣象學、天文學、海洋、航空與航天等學科和工程技術科學中發揮著越來越重要的作用,利用其可以刻畫物理科學中的基本規律,如納維-斯托克斯方程反映了粘性流體流動的基本力學規律。然而在復雜的場景下,偏微分方程的精確解析解很難求解出來。此時可以利用傳統數值求解方法對其進行近似求解。常用的數值求解方法包括有限差分法、有限體積法、有限元法等。這些數值方法將偏微分方程的求解區域離散成若干網格點,之后在離散的網格點上提供解值,從而得到方程的近似解。王敞亮[1]分析討論了計算機流體力學中控制方程有限差分法數值求解,具體從有限差分法的概念及基本思路進行了分析。廖敦明等[2]利用有限差分法對鑄件凝固過程熱應力場進行了數值求解及模擬,并且開發了基于有限差分法應力場模擬系統。文中闡述可將鑄件成形過程劃分成許多的小變形增量步,并且利用基于小變形熱彈塑性理論來對每個小變形增量步進行求解。劉鐵錘等[3]總結了國內外基于有限體積法對潰壩水流模擬的研究進展,其中介紹了應用于有限體積法的非結構網格、有限體積法中方程離散的步驟以及法向數值通量計算等。成思源[4]詳細闡述和討論了有限元法的原理及產生方法論。

雖然這些傳統數值求解方法能夠對于一些難以求出解析解的復雜的偏微分方程進行求解,但是當需要求解精度越高時,計算速度會非常緩慢,并且耗費的內存及磁盤空間等計算資源都非常驚人,差商所依賴的網格節點也會越多,此時會使得邊界條件也不好處理。

最近幾年,隨著卷積神經網絡(CNNs)的迅速發展,其在加速偏微分方程數值求解方面取得了重大突破。例如Bar-Sinai Y等[5-6]將CNN嵌入傳統數值求解過程中,用它代替有限差分系數的計算,而其余的步驟仍使用傳統方法,如計算通量、做時間步進。機器學習組件與傳統數值求解方案緊密集成,并通過端到端可微編程框架進行訓練。Kochkov D等[7]提出了一種新型的數值求解器,它不平均未解析的自由度,而是使用離散方程,在未解析網格上給出逐點精確解。Sirignano等[8]提出了一個深度學習偏微分方程模型(DPM),它是一種基于深度學習的偏微分方程增強方法,該方法優化了偏微分方程的函數形式。Um K等[9]針對偏微分方程離散化過程中出現的數值誤差,離散化誤差可以通過神經網絡進行學習,訓練好的網絡可以在本地進行評估,從而可以減少數值誤差,改進偏微分方程求解器的解。Pathak J等[10]引入了一種新的ML-PDE混合體系結構,該技術可以修正高雷諾數下湍流的粗網格模擬引起的數值誤差,同時恢復高分辨率場的估計,并且能在以低分辨率求解系統PDE的同時獲得高分辨率解軌跡。Tompson J 等[11]提出了一種數據驅動方法,即利用深度學習的近似能力和標準求解器的精度來獲得快速和高度真實的模擬。Obiols-Sales O等[12]提出了一個深度學習和物理模擬耦合的框架(CFDNet),實驗結果表明使用CFDNet顯著加快了數值求解的速度。

雖然這些方法在偏微分方程數值求解加速領域都取得了重大的突破,但是這些方法這些方法仍存在兩個關鍵的缺點:缺乏物理先驗知識的指導。多數方法建模時由于沒有考慮物理先驗知識的指導,在預測方程解的過程中忽略了部分物理信息,從而導致求解精度下降。忽略了多層級空間聯結關系和全局信息。多數方法利用神經網絡進行數值求解加速時,建模過程中采用扁平化建模,沒有進行結構化的建模,因此并沒有將多層級復雜空間聯結關系考慮在內。除了目標節點周圍節點會對該點產生影響,其他節點也會影響到目標節點,而多數方法僅僅考慮了目標節點周圍的節點,忽略了全局空間的影響,從而引起全局信息缺失。

針對以上問題,本文利用經典圖模型構建空間圖結構,所有空間網格點作為圖結點,并基于物理先驗知識確定圖結點間關系,構造融合全局信息的圖鄰接矩陣。因此,本文提出了基于圖卷積的空間離散化數值求解加速方法,并利用該方法來加速求解一維平流方程。

1 方法

1.1 一維平流方程

平流方程指的是任意變量水平輸送導致該變量的局地變化的微分方程,它是大氣運動中重要過程之一,因此,對平流方程的研究對大氣運動的研究具有非常重要的意義。平流方程最簡單的形式是一維平流方程,它被寫成:

(1)

(2)

式中:{x1,…,xN}為空間網格點;Ci+j為點xi+j處的濃度;{α-k,…,αk}為預定義的有限差分系數。

利用歐拉格式進行數值求解可以在高階情況下得到高精度的解,但是其計算速度極其緩慢,并且在復雜速度場下,如湍流速度場下,還會引起極大的數值擴散。

1.2 物理先驗知識

在進行一維平流方程空間離散化數值求解時,將方程的連續的定解區域進行離散,由于本文采用6點有限差分模板,因此之后將當前時刻每個點處的空間導數由預測得到空間離散化有限差分系數與上一時刻的該點處的濃度值及該點周圍的6個點處的濃度值相結合來表示,由此可得到當前時刻每個點處的濃度值。6點有限差分模板是借鑒的傳統數值求解方法,之后本文在此基礎上進一步根據物理場的性質定義并設計了物理先驗知識,從而對初始系數進行了設定。因此本文的物理先驗知識分為兩種:①對于每個節點周圍的6個點,離目標節點越近的點對其影響越大;②除了這6個點以外,其他節點也會對目標節點產生影響。在通過圖模型對一維平流方程空間離散化結構建模時,基于這兩個物理先驗知識確定圖結點間關系,從而構造融合全局信息的圖鄰接矩陣。

基于這兩個物理先驗知識,通過圖模型構建空間結構鄰接矩陣,根據先驗知識①,對于每個節點周圍的6個點,在鄰接矩陣上將離目標節點更近的點的值設置的更大;根據先驗知識②,由于除了每個點周圍的6個點外的其他點也會對目標節點產生影響,因此不應將其他節點一概設置為0,將它們設置為與0接近的值。例如下面的矩陣。

(3)

1.3 基于圖卷積的一維平流方程空間離散化數值求解加速方法

(前向求解過程中,使用圖卷積物理神經網絡(紅色框部分)預測出每一層的有限差分系數,之后去計算空間差分、時間步進,從而得到每一層的濃度值。在訓練時,使用由基線求解器求解得到的高分辨率解下采樣后的低分辨率數據集對模型進行優化。 In the forward solution process, the physical neural network of graph convolution (red box part) is used to predict the finite difference coefficient of each layer, and then the spatial difference and time step are calculated to obtain the concentration value of each layer. During training, the low resolution data set sampled from the high resolution solution obtained by the baseline solver is used to optimize the model.)圖1 基于圖卷積的一維平流方程空間離散化數值求解器整體框架Fig.1 Overall framework of spatial discretization numerical solver for one-dimensional advection equation based on graph convolution

1.3.2 GCPNN部分 本文利用GCPNN(圖1紅色框部分)來預測出有限差分系數,在該部分使用由第一層濃度值得到的輸入特征矩陣(Feature matrix)H(1)和由網格節點間的空間結構關系得到的鄰接矩陣(Adjacency matrix)A相乘,之后將得到的結果經過全連接層,從而得到第二層的特征,接著再將該層的特征與鄰接矩陣相乘后經過全連接層,得到第三層的特征,重復上述操作,直到得到第n層的特征,之后再經過CNN卷積操作。

在最后加入了多項式精度約束(Polynomial constraint),它可以強制機器學習預測系數滿足n階多項式約束,從而使近似誤差衰減為O(Δxn),如果學習的離散化適合于網格尺度上平滑的解,多項式精度約束將確保可以恢復經典的有限差分方法。

1.3.3 構建圖結構及鄰接矩陣 基于圖模型的一維平流方程空間離散化結構建模需要構建網格節點間的圖結構關系,首先確定圖結構的結點和邊。

結點:將空間離散化后的高分辨率網格點下采樣后得到的低分辨率點作為圖結構的結點。

邊:基于物理先驗知識通過距離遠近空間多層級化來設計圖結構的邊的權重。

之后建立所有網格點的空間結構關系;具體如下:假設空間離散化后共有m個網格點,每個網格點都有濃度值,每個網格點的濃度值由上一層該點處的值和該點周圍2k個網格點的值來表示(即該點左邊和右邊各k個點。這里k是經實驗得到的最優周圍點數量,即通過k個周圍點預測得到的中心點濃度值與真解間誤差最小,在本文中取k為3,即6點有限差分模板。)。根據每個網格點的濃度值得到輸入特征矩陣H(1),H(1)形狀為(m,1),如下為特征矩陣H(1),其中C1,…,Cm分別為m個網格點上的濃度值。

(4)

由空間結構關系可得到鄰接矩陣A(m×m)為:

(5)

(6)

其中,H為預測出的每一時間層的特征矩陣,輸入層為H(1);σ是非線性激活函數。

1.3.4 多步損失函數 在訓練過程中,本文使用多步損失函數,相對于單個時間步,多步損失函數可以穩定前向積分。多步損失函數即累積預測10個時間步的前向求解器預測結果,然后在整個空間域上與該時間段內的真解計算平均絕對誤差(MAE)來做為損失函數:

(7)

其中:N表示10個時間步;M表示每個時間步上的網格節點個數。

2 實驗過程及實驗結果

2.1 數據集

本文使用帶有單調通量限制器的二階VanLeer平流方案來作為基線求解器,之后以64×64的分辨率來運行基線求解器來獲得參考真解(ground truth)。之后對獲得的較高分辨率結果進行下采樣得到32×32低分辨率的結果,即可為以粗糙網格訓練的機器學習模型生成訓練和測試數據集。

2.2 訓練優化方法

本文使用Tensorflow框架實現模型的構建,并采用8 000個解作為訓練樣本,2 000個解作為測試樣本。在訓練時,迭代次數(epoch)設置為300,批量大小(batch size)設置為128。

2.3 實驗指標

(8)

其中N表示網格節點的個數。

2.4 實驗設計

2.4.1 對比實驗 本文將基于圖卷積模型的求解器與基線求解器、一階Vanleer平流格式求解器、基于cnn模型的求解器進行對比,將預測出的解與真解(Ground truth)之間的平均絕對誤差(MAE)作為對比指標,圖2展示了在這些求解器上所有測試樣本的平均絕對誤差(MAE)隨時間的變化,圖卷積模型的誤差最小且比基線二階VanLeer格式求解器小8倍,并且圖卷積模型比cnn模型的誤差小。

圖2 在測試樣本上得到的誤差Fig. 2 Error obtained on the test sample

在測試樣本上圖卷積模型單步預測時間為0.35 s,cnn模型單步預測時間為0.29 s。在經過了250個時間步后,使用CNN模型的平均絕對誤差的值約為0.007,使用圖卷積模型的平均絕對誤差的值約為0.002。因此雖然圖卷積單步預測時間與cnn模型相比多消耗了20%,但精度卻比CNN模型提高了71%。

2.4.2 消融實驗 (1)基于物理先驗知識①的鄰接矩陣值更換

為了驗證離某個網格節點越近的點對其影響越大的假設成立,本文在之前的鄰接矩陣基礎上進行了改動,且設計出了三種鄰接矩陣,即將矩陣中的6個1替換為另外6個接近于1的數,且越靠近該節點的點處的值越大,三種模板設計分別為[0.97,0.98,0.99,0.99,0.98,0.97]、[0.75,0.85,0.95,0.95,0.85,0.75]、[0.7,0.8,0.9,0.9,0.8,0.7],并用這三種模板得到的鄰接矩陣去進行訓練,訓練得到的test loss如表1所示。從表中可以看出使用這三種鄰接矩陣訓練得到的test loss與使用原鄰接矩陣訓練得到的test loss相差無幾,且使用[0.75,0.85,0.95,0.95,0.85,0.75]這種模板得到的test loss最低,由此可以得出離某個網格節點越近的點對其產生的影響越大,且鄰接矩陣使用[0.75,0.85,0.95,0.95,0.85,0.75]這種模板來進行設計所得結果會更好。

表1 不同鄰接矩陣模板訓練結果對比Table 1 Comparison of training results of different adjacency matrix templates

(2)基于物理先驗知識②的鄰接矩陣值更換

為了驗證對于某個網格節點,不僅該點周圍的6個點對其有影響,其他點也會對該點產生影響,本文又設計了四種鄰接矩陣,這次替換原鄰接矩陣中的0,分別為①將所有0替換為0.1;②將所有0替換為0.2;③將所有0替換為0.001;④將所有0替換為0.003。并使用這四種鄰接矩陣去進行訓練,訓練得到的test loss如表2所示。從表中可以看出使用將所有0替換為0.001得到的test loss最低,由此可以看出除了某個點周圍的6個點對該點有影響,其他點也對該點會產生影響,且鄰接矩陣使用將所有0替換為0.001來進行設計所得結果會更好。

表2 不同鄰接矩陣模板訓練結果對比Table 2 Comparison of training results of different adjacency matrix templates

(3)有限差分模板更改實驗

為了得到更好的有限差分模板,本文又設計了三種鄰接矩陣,這次改變原鄰接矩陣中每行1的個數(原鄰接矩陣中為6),分別將其改為 5、7、8個,使用這些鄰接矩陣訓練得到的test loss如表3所示,從表中可以看出原鄰接矩陣比其他三種矩陣訓練得到的test loss都低,由此可以看出,6點有限差分模板是更好的。

表3 不同鄰接矩陣模板訓練結果對比Table 3 Comparison of training results of different adjacency matrix templates

3 結論與展望

本文針對一維平流方程,在目前現有數據驅動離散化方法的基礎上創新性的提出了基于圖卷積的一維平流方程空間離散化數值求解加速方法,該方法在物理先驗指導下基于圖模型對一維平流方程空間離散化進行結構建模,并利用圖卷積算法構建出預測模型,利用該模型學習最優的有限差分系數,實現加速數值求解。之后并進行了實驗分析,包括對比實驗和消融實驗,主要結論如下:

(1)本文提出的基于圖卷積的一維平流方程空間離散化數值求解加速方法不僅可以在低分辨率下得到高精度的解,并且能加速求解,同時該方法的預測結果顯著優于基線求解器和CNN模型求解器。

(2)根據物理先驗知識①進行鄰接矩陣值更換,可得出結論即離某個網格節點越近的點對其產生的影響越大,且鄰接矩陣使用[0.75,0.85,0.95,0.95,0.85,0.75]這種模板來進行設計所得結果會更好。

(3)根據物理先驗知識②進行鄰接矩陣值更換,可得出結論即除了某個點周圍的6個點對該點有影響,其他點也對該點會產生影響,且鄰接矩陣使用將所有0替換為0.001來進行設計所得結果會更好。

(4)通過更改有限差分模板中1的個數來進行實驗可得出結論即使用6點有限差分模板達到的效果最好。

后續工作將考慮將本文方法應用到其他數值離散方案上,如有限體積法、有限元法、偽譜方法等。另外,在后續工作中還將考慮將本文方法推廣到高維平流方程數值求解問題。